bugGNU Octave - Bugs: bug #38143, lsode gives incorrect results for...

 
 

bug #38143: lsode gives incorrect results for basic linear vibration problems

Submitter:  Mario <nopria>
Submitted:  Tue 22 Jan 2013 04:20:02 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Invalid / Not an Octave Bug Assigned to:  None
Originator Name:  Mario Open/Closed:  * Closed
Release:  * 3.6.1 Operating System:  * Microsoft Windows
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Mon 04 Mar 2013 07:15:12 PM UTC, comment #2: 

Please close this report, thanks!

Juan Pablo Carbajal <juanpi>
Group Member
Tue 22 Jan 2013 04:46:23 PM UTC, comment #1: 

When I run your function and plot y(:,1) vs t, I get something that looks like your maxima solution.

John W. Eaton <jwe>
Group administrator
Tue 22 Jan 2013 04:20:02 PM UTC, original submission:  

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));


Mario <nopria>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #27311:  lsode solution.jpg added by nopria (63KiB - image/jpeg)
file #27310:  analytic.jpg added by nopria (44KiB - image/jpeg)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jordigh (Updated the item)
  • -email is unavailable- added by juanpi (Posted a comment)
  • -email is unavailable- added by rik5 (Updated the item)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by nopria (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2013-03-04 jordigh StatusWorks For Me Invalid / Not an Octave Bug
        Open/ClosedOpen Closed
    2013-02-01 rik5 StatusNone Works For Me
    2013-01-22 nopria Attached File- Added analytic.jpg, #27310
        Attached File- Added lsode solution.jpg, #27311
        Attached File- Added manual Adam-Moulton implementation.jpg, #27312

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code