Solving ODE numerically in Mathematica - I get 'ndnum' error

  • #1

Main Question or Discussion Point

Hello everyone!
I'm trying hard to solve numerically a system of coupled differential equations of first order, but I get this error everytime.. I can't find the reason.. maybe you can help me, I'd really apreciate that.

This is the code:

Jg=0.000043;
Kg=0.5;
Bg=0.06;
Bpm=0.5;
r=0.11;
M=0.0000677;
N1=100;
N2=-200;
t1a=9.7;
t1d=0.2;
T1=20;
t2a=2.4;
t2d=1.9;
T2=26;
kte=0.4;
ktl=0.76;
ks=1.2;
ftc=1.5;
F0=0.7;
l0=0.2;
c=2;
kml=0.7;
kme=1;
kpm=0.32;
lmc=1.2;
lms=1.4;

Fact[x_,y_]=F0*(1-(((x-1)/0.5)^2))*((l0-y)/(l0+c*y));
Kt[z_]=(kte*z)+ktl;
Fpe[x_]=(kml/kme)*(exp ((kml*(x-lms)*(180/\[Pi]*r))-1));

S= NDSolve[{x1'[t]==x2[t],x2'[t]==(1/Jg)*(x7[t]-x8[t]-(Bg*x2[t])-(Kg*x1[t])),x3'[t]==x4[t],x4'[t]==(980/M)*(x7[t]-(x9[t]*Fact[x3[t],x4[t]])-Fpe[x3[t]]-(Bpm*(180/(N[\[Pi]]*r))*x4[t])),x5'[t]==x4[t],x6'[t]==(980/M)*(x8[t]-(x10[t]*Fact[x5[t],x6[t]])-Fpe[x5[t]]-(Bpm*(180/(N[\[Pi]]*r))*x6[t])),x7'[t]==Kt[x7[t]]*(-x2[t]-((180/(N[\[Pi]]*r))*x4[t])),x8'[t]==Kt[x8[t]]*(-x2[t]-((180/(N[\[Pi]]*r))*x6[t])),x9'[t]==(1/(t1a*(N[HeavisideTheta[t]-HeavisideTheta[t-T1]])+t1d*(N[HeavisideTheta[t-T1]])))*(N1-x9[t]),x10'[t]==(1/(t2a*(N[HeavisideTheta[t]-HeavisideTheta[t-T2]])+t2d*(N[HeavisideTheta[t-T2]])))*(N2-x10[t]),x1[0]==0,x2[0]==0,x3[0]==4,x4[0]==0,x5[0]==4,x6[0]==0,x7[0]==20,x8[0]==20,x9[0]==0.17,x10[0]==0.17},{x1[t],x2[t],x3[t],x4[t],x5[t],x6[t],x7[t],x8[t],x9[t],x10[t]},{t,0,5}]

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>



Every constant is numerically defined above.. What am I doing wrong here??

Thanks a lot!
 

Answers and Replies

  • #2
29,738
6,071
There may be more errors, but there are at least two. First, in your definition of Fpe you are using the undefined symbol "exp". I assume that you meant "Exp[]" (note the square brackets as well as the capitalization). Second, HeavisideTheta[t] is undefined for t=0. You may want to use UnitStep[t] instead.

There may be other problems, but frankly this is very difficult to read. You need to really test your inputs to a function like this and make sure that all of them are correct. Chances are pretty high that your constants like M and Jg will cause numerical instabilities.
 
  • #3
Thanks very much! I still got an error, but it's different now... I fixed the functions errors and adjust a little bit the constants M and Jg becouse they we're close to zero..


[tex]\text{Fact}[\text{x$\_$},\text{y$\_$}]=\text{F0}*(1-(((x-1)/0.5){}^{\wedge}2))*((\text{l0}-y)/(\text{l0}+c*y));[/tex]
[tex]\text{Kt}[\text{z$\_$}]=(\text{kte}*z)+\text{ktl};[/tex]
[tex]\text{Fpe}[\text{x$\_$}]=(\text{kml}/\text{kme})*(\text{Exp}[\text{kml}*(x-\text{lms})*(180/\pi *r)]-1);[/tex]

S= NDSolve[{
x1'[t]==x2[t],
x2'[t]==(1/Jg)*(x7[t]-x8[t]-(Bg*x2[t])-(Kg*x1[t])),
x3'[t]==x4[t],x4'[t]==(980/M)*(x7[t]-(x9[t]*Fact[x3[t],x4[t]])-Fpe[x3[t]]-(Bpm*(180/(Pi*r))*x4[t])),
x5'[t]==x4[t],
x6'[t]==(980/M)*(x8[t]-(x10[t]*Fact[x5[t],x6[t]])-Fpe[x5[t]]-(Bpm*(180/(Pi*r))*x6[t])),x7'[t]==Kt[x7[t]]*(-x2[t]-((180/(Pi*r))*x4[t])),
x8'[t]==Kt[x8[t]]*(-x2[t]-((180/(Pi*r))*x6[t])),
x9'[t]==(1/(t1a*(UnitStep[t]-UnitStep[t-T1])+t1d*(UnitStep[t-T1])))*(N1-x9[t]),
x10'[t]==(1/(t2a*(UnitStep[t]-UnitStep[t-T2])+t2d*(UnitStep[t-T2])))*(N2-x10[t]),
x1[0]==0,x2[0]==0,x3[0]==4,x4[0]==0,x5[0]==4,x6[0]==0,x7[0]==20,x8[0]==20,x9[0]==0.17,x10[0]==0.17},
{x1[t],x2[t],x3[t],x4[t],x5[t],x6[t],x7[t],x8[t],x9[t],x10[t]},{t,1,5}]


NDSolve::ndsz: At t == 0.02201988618427872`, step size is effectively zero; singularity or stiff system suspected. >>
NDSolve::ndsz: At t == 0.02201988618427872`, step size is effectively zero; singularity or stiff system suspected. >>
{}

I have separated the part of NDSolve into sentences, so it's a little bit more "readable"
What does that means now ??
Thank You!
 
Last edited:
  • #4
29,738
6,071
That means that your system is running into numerical instability problems.

First, you need to plot all of your x1', x2', etc. values across the whole range of possible x and t values to see if any of them are blowing up anywhere. This is what I meant above by testing all of your inputs.

If they are all well behaved then you may want to try the "StiffnessSwitching" controller for the Method option. You probably should also specify a very small StartingStepSize and maybe even consider increasing your WorkingPrecision.

It may be that, even with all of the above, your system is simply a ill-conditioned system.
 
  • Like
Likes Dr.ahmad adnan

Related Threads on Solving ODE numerically in Mathematica - I get 'ndnum' error

Replies
4
Views
5K
Replies
0
Views
2K
Replies
3
Views
3K
  • Last Post
Replies
6
Views
3K
Replies
4
Views
6K
  • Last Post
2
Replies
34
Views
2K
  • Last Post
Replies
4
Views
2K
Top