help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

lsode precision


From: Mario Maio
Subject: lsode precision
Date: Mon, 21 Jan 2013 17:03:08 +0100
User-agent: Mozilla/5.0 (Windows NT 6.1; WOW64; rv:17.0) Gecko/20130107 Thunderbird/17.0.2

I use lsode to solve a simple standard vibration problem with Octave 3.6.1:

function ydot=f(y,t) % 2nd order equation conversion to a 1st order system
    ydot(1)=y(2);
    ydot(2)=eqd(t,y(1),y(2));
end

function xdd = eqd(t,x,xd) % generic vibration
    m=1;
    K=1;
    omega=sqrt(K/m); % pulsazione naturale oscillatore
    fattore_smorzamento=0.07;
    S=2*fattore_smorzamento*m*omega;
    ampiezza_forzante=1;
    Omega=2.5; % pulsazione forzante
    forzante=ampiezza_forzante/m*sin(Omega*t);
    xdd=(-S*xd-K*x+forzante)/m;
end

y=lsode("f",[0;0],(t=(0:0.1:20*pi)'));
figure(1);plot(t,y(:,2))


But the result is qualitatively different from a manual Adam-Moulton corrector implementation which agrees with the text book the example was taken from.

Add the following code to the previous:

function [x,xd,xdd]=pred_corr(dt,tp,xp,xdp,xddp,tol)
    % tp, xdp, xddp are values of t, xd and xdd at previous step
    xddg=xddp; % as initial guess xdd is taken same as xddp
    for j=1:100
        xd=xdp+dt/2*(xddp+xddg);
        x=xp+dt/2*(xdp+xd);
        xdd=eqd(tp+dt,x,xd);
        if abs(xdd-xddg)<tol
            break;
        else
            xddg=xdd;
        end
    end
end

dt=.1;
tv=0:dt:20*pi;
t=tv(1);

x=0; % initial values
xd=0;

xdd=eqd(t,x,xd);

xv=[x];
for t=tv(1:end-1)
    [x,xd,xdd]=pred_corr(dt,t,x,xd,xdd,.000001);
    %display([t+dt x xd xdd])
    xv(end+1)=x;
end

figure(2);plot(tv,xv)


Is lsode unprecise? I have to change some settings to get more precision from lsode? (actually I already tried to change some settings but the result is the same)

Thanks.


reply via email to

[Prev in Thread] Current Thread [Next in Thread]