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

[LeoDimieri]

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!

 

[Dale]

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.

 

[LeoDimieri]

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..\text{Fact}[\text{x$\_$},\text{y$\_$}]=\text{F0}*(1-(((x-1)/0.5){}^{\wedge}2))*((\text{l0}-y)/(\text{l0}+c*y));
\text{Kt}[\text{z$\_$}]=(\text{kte}*z)+\text{ktl};
\text{Fpe}[\text{x$\_$}]=(\text{kml}/\text{kme})*(\text{Exp}[\text{kml}*(x-\text{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/(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!

 

[Dale]

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.