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

 
 

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

Submitted by:  Mario <nopria>
Submitted on:  Tue 22 Jan 2013 04:20:02 PM UTC  
 
Category: LibrariesSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: InvalidAssigned to: None
Originator Name: MarioOpen/Closed: Closed
Release: 3.6.1Operating System: Microsoft Windows

Add a New Comment (Rich MarkupRich Markup):
   

You are not logged in

Please log in, so followups can be emailed to you.

 

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

Please close this report, thanks!

Juan Pablo Carbajal <juanpi>
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>
Project 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=2fattore_smorzamentom*omega;
ampiezza_forzante=1;
Omega=2.5; % pulsazione forzante
forzante=ampiezza_forzante/msin(Omegat);
xdd=(-Sxd-Kx+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^(-omegatxi)(((4omega^2Omegaxi^2+2Omega^3-2omega^2Omega)sin((tsqrt(4omega^2-4omega^2xi^2))/2)F)/(4momega^2Omega^2xi^2sqrt(4omega^2-4omega^2xi^2)+mOmega^4sqrt(4omega^2-4omega^2xi^2)-2momega^2Omega^2sqrt(4omega^2-4omega^2xi^2)+momega^4sqrt(4omega^2-4omega^2xi^2))+(2omegaOmegaxicos((tsqrt(4omega^2-4omega^2xi^2))/2)F)/(4momega^2Omega^2xi^2+mOmega^4-2momega^2Omega^2+momega^4))-
((2omegaOmegacos(Omegat)xi+(Omega^2-omega^2)sin(Omegat))F)/(4momega^2Omega^2xi^2+mOmega^4-2momega^2Omega^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 File(s):
   
   
Comment:
   

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

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 6 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 04 Mar 2013 07:34:43 PM UTCjordighStatusWorks For Me=>Invalid
      Open/ClosedOpen=>Closed
    Fri 01 Feb 2013 12:01:15 AM UTCrik5StatusNone=>Works For Me
    Tue 22 Jan 2013 04:20:02 PM UTCnopriaAttached 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.1-cleanup