Re: lsode problem with simple ODE system.

[Olaf Till]

On Sat, Aug 31, 2013 at 12:05:33PM -0700, Paolo Ferraresi wrote:
> Hello, I have a problem with "lsode".
> The following script, is solved in seconds with a simple fixed step
> Runge-Kutta (23 order) while "lsode" seems not able to integrate this ode
> system.

With fixed step you get some result for anything, as long as the user
function returns some real values for xdot. But the result is not
necessarily meaningful.

With stepsize controled by a guess of accuracy, as in lsodes
algorithm, there may be problems if something in the user function
xdot is just wrong, or if there are discontinuities in its returned
value.

In your case, the use of 'sign()' (twice in your function) is likely
to cause discontinuities. And if you model some real world problem,
I'd suspect that the use of 'sign()' is not correct, too. Try to
replace it with some general term that "captures the change in sign
internally" in some way ... (just leaving out the two 'sign()' calls
makes integration successful for your case, but I don't know how
meaningful xdot is then).

Sometimes lsode just needs a long time, but not typically with only
three differential equations. If lsode gets stuck at a particular
value of 't', typically there is some mistake in the differential
equations.

Olaf

> %% temp.m
> 1;
> 
> global beta g m mu mu_d V Ap kd km kr F0 xmax Al Cr;
> g = 9.81;
> m = 0.1;
> mu_d = 0.1; 
> Dp = 24e-3
> Ap = pi*Dp^2/4;
> km = 40.151e+3
> kr = 1e+6;
> L0 = 60.9e-3
> L1 = 47e-3
> L2 = 45e-3
> xmax = 7e-3;
> F0 = -km*(L0-L1)
> 
> Cd = 0.7;
> rho = 920;
> ds = 0.8e-3;
> As = pi*ds^2/4;
> kd = Cd*As*sqrt(2/rho);
> beta = 1.75e+9;
> V = 1e-3;
> 
> L = 5e-3;
> Al = pi*Dp*L;
> Cr = 0.02e-3;
> mu = 0.039;
> 
> %% z(1) = x, z(2) = v, z(3) = p
> function qdot = q(z,t)
>       global beta g m mu mu_d V Ap kd km kr F0 xmax Al Cr;
>       x = z(1);
>       v = z(2);
>       p = z(3);
>       
>       pdot = (beta/V)*(-kd*sqrt(abs(p))*sign(p)-Ap*v);
>       Fm = F0-km*x;
>       Fc = 0.0;
>       if (x<0)
>               Fc = -x*kr;
>               elseif (x>xmax)
>                       Fc = (xmax-x)*kr;
>               endif
>       Fr = -mu_d*m*g*sign(v); %% attrito radente
>       Fv = -mu*Al*v/Cr;
>       F = p*Ap+Fm+Fc+Fr+Fv;
>       vdot = F/m;
>       xdot = v;
>       
>       qdot = [ xdot ; vdot ; pdot ];
>       endfunction
> 
> %t = linspace(0,1,1000000);
> %x = lsode("q", [ xmax ; 0 ; 0 ] , t);
> [t,x] = rk2fixed(@q,[0 0.5],[xmax ; 0 ; 0], 10000 , 1 , 0 , 0);
> 
> 
> 
> 
> --
> View this message in context: 
> http://octave.1599824.n4.nabble.com/lsode-problem-with-simple-ODE-system-tp4656959.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://mailman.cae.wisc.edu/listinfo/help-octave

-- 
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc
Description: Digital signature