lsode problem with simple ODE system.

[Paolo Ferraresi]

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.
Even if you change the parameters ("stiff", "adams", working with tolerances
and other parameters). With "stiff" i got an error. With "adams" it takes an
infinite time ...

Who can help me. I would prefer to use the native "lsode" rather than rk.

--------------------------------------------------------------------------------

%% 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.