Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Mar 17, 2015 #1
    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!
     
  2. jcsd
  3. Mar 17, 2015 #2

    Dale

    Staff: Mentor

    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.
     
  4. Mar 17, 2015 #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: Mar 17, 2015
  5. Mar 17, 2015 #4

    Dale

    Staff: Mentor

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Solving ODE numerically in Mathematica - I get 'ndnum' error
Loading...