Help weird NDSolve error, system of ODE and an algebraic equation

Questions
, , , ,
1

Hi, I have a system of 2 ODEs and 2 algebraic equations. I tried to set the 2 algebraic equations as ODEs but I’m not sure thats completely correct. This system of equations is supposed to model a biological system, and the 2 algebraic equations start at some equilibrium-value that I don’t really know, so I want mathematica to calculate it using those equations. Theoretically, since the system will always be at equilib, it won’t shift, so I dont think i can use them as differential equations. Anyways this is the error:

NDSolve has computed initial values that give a zero residual for the
\ differential-algebraic system, but some components are different
from \ those specified. If you need them to be satisfied, giving
initial \ conditions for all dependent variables and their derivatives
is \ recommended.

I don’t get it because I (think) I have defined all initial conditions. maybe it doesn’t like the algebraic equations?
I have attached the file. I am using mathematica 9 on win8.
Any help would be awesome, thanks!!

2

Your deq4 looks like deq4 := r[t] == -p[7] r[t]; This says that -p[7] is one? But you later on write p[7] -> 0.003

It will be better to post the actual code here than attaching a notebook. The code is not large.

I also do not know why you used indexed variables p[3] etc.. instead of normal variables p3 etc… There is no need at all to use indexed variable instead of normal variables. Having [] all over the code makes it hard to read.

I would also like to suggest not to use l for variable name, as it looks like a 1 on many screen and fonts.
The error you get indicates an inconsistent initial condition. Please see NDSolveIntroductoryTutorialDAEs

3

Here is my new code:

    DAE = (s'[
         t] == -p3 (r[t] s[t] c[t]/(p1*p2))/((1 + s[t]/p1) (1 + c[t]/p2))
         w'[t] 
w'[t] == 
        p3*p5 (r[t] s[t] c[t]/(p2*p2))/((1 + s[t]/p2) (1 + c[t]/p3)) - w[t]
         
 c[t] == -(1 + p1/s[t]) (r[t]*s[t]*
             c[t]/(p1*p2))/((1 + s[t]/p1) (1 + c[t]/p2)) - p7*c[t]
         
 r[t] == -p7*r[t]);

sol = NDSolve[{DAE, s[0] == 150, w[0] == 0} /. {p1 -> 50, p2 -> 70.5, 
    p3 -> .16, p5 -> 1600, p[6] -> .05, p7 -> 0.03}, {s, w, c, r}, {t,
    0, 120}]

And now I get this error:

NDSolve::indexss: The DAE solver failed at t = 0.0016`. The solver is
intended for index 1 DAE systems and structural analysis indicates
that the DAE is structurally singular.

Now I really have no idea what it means. Help please?

4

There are a couple syntax errors in your code. I would not recommend assigning a set of equations to a variable and then passing that variable into an NDSolve or DSolve. Additionally the constants that are present in your differential equations do not have a common notation – i.e. p[6] instead of p6. And I cannot see that p6 is called in your set of equations. Is that meant to be the case?
From a theoretical perspective, you should be looking for analytic solutions before brute numerical solutions. Have you tried that?

A cleaner version of the code to search for analytic solutions might look like:

      p1 = 50;
      p2 = 70.5;
      p3 = .16;
      p5 = 1600;
      p6 = .05;
      p7 = 0.03;
      sol1 = DSolve[{s'[t] == -p3 (r[t] s[t] c[t]/(p1*p2))/((1 + s[t]/p1) (1 + c[t]/p2)) , 
      w'[t] ==p3*p5 (r[t] s[t] c[t]/(p2*p2))/((1 + s[t]/p2) (1 +c[t]/p3)) - 
      w[t], c[t] == -(1 + p1/s[t]) (r[t]*s[t]*
     c[t]/(p1*p2))/((1 + s[t]/p1) (1 + c[t]/p2)) - 
 p7*c[t] r[t] == -p7*r[t], s[0] == 150, w[0] == 0, c[0] == 0.05,

r[0] == 0.05}, {s[t], w[t], c[t], r[t]}, {t}]

To look for numerical solutions, a cleaner version of the last section of the code might look like:

sol = NDSolve[{s'[
  t] == -p3 (r[t] s[
      t] c[t]/(p1*p2))/((1 + s[t]/p1) (1 + c[t]/p2)) w'[t], 
w'[t] == 
 p3*p5 (r[t] s[t] c[t]/(p2*p2))/((1 + s[t]/p2) (1 + c[t]/p3)) - 
  w[t], c[t] == -(1 + p1/s[t]) (r[t]*s[t]*
      c[t]/(p1*p2))/((1 + s[t]/p1) (1 + c[t]/p2)) - 
  p7*c[t] r[t] == -p7*r[t], s[0] == 150, w[0] == 0} /. {p1 -> 50, 
p2 -> 70.5, p3 -> .16, p5 -> 1600, p[6] -> .05, p7 -> 0.03}, {s, 
 w, c, r}, {t, 0, 120}]