Solving simultaneously algebraic and differential equations

Hello,

I'm modelling a polymerization process and need to solve both kind of equations at the same time.

I have created until now, 2 .m files:

1)contain the diferential equations that need to be solved

2)contain 1 implicit equation to find his root

The equations to represent the system are like this:

Mw=f(x,Vp,Nd,Mw)
Vp'=f(N1,Mw)
Nd=f(N,Vp)
N1=f(N,N0)
N0'=f(N0,N1)
N'=f(Mw,Nd)
x'=f(N1,Mw)
...

For example to calculate dVp/dt i need Mw, and to calculate Mw i need Vp, so both equations need to be solved simultaneously.

The structure of the file to solve the diferential equations is:

function [f]=sys(t,y)

function [f]=sys(t,y)
IntMch=... %Here i define the parameters
Vp=0; DVp=0.00001; y(1)=1;
While Vp<y(1)
Vp=Vp+DVp;
Mw = fzero(@(Mw)concmw(Mw,Vp,Nd,...),10^-3);
f(1) = Mw*...+...;
f(2) = y(2)*...+...;
f(3) = y(2)*...+y(3)*...
f=f';
end;

The file to solve concmw is:

function F=concmw(Mw,Vp,Nd,...)
F=(Mw/...)^(1/0.6)*Vp+Mw*Nd*...+Mw*...

Finally, on the command window, i call the first file with:

tr=[0,7200];x0=[10^-5,0,0];[t,x]=ode45('sys',tr,x0)

When i try to run, i get this:

Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at -0.00028 is -0.034758-0.0012509i.)
Check function or try again with a different starting value.

The problem is the

^(1/0.6)

But i get this message on each iteration and never stops.

So, there is a way to solve this kind of system? How to restrict the value of Mw to avoid that function goes to a complex value?

Later like in Mw=f(Vp,Mw,Nd,...) I need to do the same procedure to get Nd=f(N,Vp,...), so all the system must get a valid result.

Sorry if i writed wrong sentences, english is not my native language.

Thanks in advance.

Answer 1

[Walter Roberson]

In your subexpression (Mw/...)^(1/0.6) if the values not shown are positive, then when Mw becomes negative, you are going to be taking the root of a negative number, which is going to give you a complex value.

Mw is your argument to fzero(), and in the form if fzero() you are using, there is no way to prevent fzero() from trying negative Mw.

You could try using abs(Mw) . That would present a problem in theory as abs(Mw) is not a continuous function. sqrt(Mw^2) would provide the necessary theoretical continuity, but since it is doing numeric work anyhow rather than symbolic, might as well use abs().

The alternative is to use the alternate form of fzero() in which your provide x0 as a pair of values instead of a single value. The pair provide upper and lower bounds on the search. It is important to fzero() that the signs of the function be different at the two endpoints.