[Octave-bug-tracker] [bug #38143] lsode gives incorrect results for basic linear vibration problems

[Mario]

URL:
  <http://savannah.gnu.org/bugs/?38143>

                 Summary: lsode gives incorrect results for basic linear
vibration problems
                 Project: GNU Octave
            Submitted by: nopria
            Submitted on: mar 22 gen 2013 16:20:02 GMT
                Category: Libraries
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Incorrect Result
                  Status: None
             Assigned to: None
         Originator Name: Mario
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 3.6.1
        Operating System: Microsoft Windows

    _______________________________________________________

Details:

I spent an hour writing a nice and complete bug report, but everything was
deleted when I clicked the link to answer the final question!!
Sorry but I won't write again that nice report, I will just forward my
original post on the forum to let you know that lsode gives incorrect results
for basic linear vibration problems.

********************************************************************************
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)
*****************************************************************************
Here is the analytic solution found with Maxima (copy and paste it):

x:%e^(-omega*t*xi)*(((4*omega^2*Omega*xi^2+2*Omega^3-2*omega^2*Omega)*sin((t*sqrt(4*omega^2-4*omega^2*xi^2))/2)*F)/(4*m*omega^2*Omega^2*xi^2*sqrt(4*omega^2-4*omega^2*xi^2)+m*Omega^4*sqrt(4*omega^2-4*omega^2*xi^2)-2*m*omega^2*Omega^2*sqrt(4*omega^2-4*omega^2*xi^2)+m*omega^4*sqrt(4*omega^2-4*omega^2*xi^2))+(2*omega*Omega*xi*cos((t*sqrt(4*omega^2-4*omega^2*xi^2))/2)*F)/(4*m*omega^2*Omega^2*xi^2+m*Omega^4-2*m*omega^2*Omega^2+m*omega^4))-
((2*omega*Omega*cos(Omega*t)*xi+(Omega^2-omega^2)*sin(Omega*t))*F)/(4*m*omega^2*Omega^2*xi^2+m*Omega^4-2*m*omega^2*Omega^2+m*omega^4);
m:1;
omega:1;
xi:0.07;
Omega:2.5;
F:1;
load("draw");
draw2d(explicit(x,t,0,50));





    _______________________________________________________

File Attachments:


-------------------------------------------------------
Date: mar 22 gen 2013 16:20:02 GMT  Name: analytic.jpg  Size: 44kB   By:
nopria

<http://savannah.gnu.org/bugs/download.php?file_id=27310>
-------------------------------------------------------
Date: mar 22 gen 2013 16:20:02 GMT  Name: lsode solution.jpg  Size: 63kB   By:
nopria

<http://savannah.gnu.org/bugs/download.php?file_id=27311>
-------------------------------------------------------
Date: mar 22 gen 2013 16:20:02 GMT  Name: manual Adam-Moulton
implementation.jpg  Size: 52kB   By: nopria

<http://savannah.gnu.org/bugs/download.php?file_id=27312>

    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?38143>

_______________________________________________
  Messaggio inviato con/da Savannah
  http://savannah.gnu.org/