Re: lsode precision

[Juan Pablo Carbajal]

Mario,

It seems you are solving a linear problem. Can you integrate the
equations analytically and compare with a exact analytic solution?
That is indeed one way of testing numeric integrators, you compare
their solutions to know exact solutions.

It seems your equations are

d²x/dt² + (S/m) * dx/dt + (K/m) * x = (A/m) * sin (Omega*t)

If so, you can find the steady state analytic solution here
http://en.wikipedia.org/wiki/Harmonic_oscillator#Sinusoidal_driving_force

Using the equations there and your data I get for the amplitude of the
steady state
m     = 1;        % Mass
w0   = 1;        % Natural frequency
zeta = 0.07;  % Damping ratio
w     = 2.5;     % Forcing frequency
F0    = 1;        % Forcing amplitude
Zm = sqrt ( (2*w0*zeta)^2 + (w0^2 - w^2)^2 / w^2 );

A = F0 / (m*Zm*w);

A =  0.19005

and from the last 20 seconds of your simulations (using your manual
implementation of the solver) I get
A = 0.189

Also integrating with ide45 from odepkg
[t_ode y_ode]=ode45(@(t,x)f(x,t),[0 20*pi],[0;0]);

Gives the same solution as your manual implementation

So it seems that lsode is not solving the problem correctly. Could you
report a bug and see what other people say?
https://savannah.gnu.org/bugs/?group=octave



On Tue, Jan 22, 2013 at 6:03 AM, Ed Meyer <address@hidden> wrote:
>
>
> On Mon, Jan 21, 2013 at 8:03 AM, Mario Maio <address@hidden> wrote:
>>
>> 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.
>> _______________________________________________
>> Help-octave mailing list
>> address@hidden
>> https://mailman.cae.wisc.edu/listinfo/help-octave
>
>
> Mario,
> I ran your job and it looks plausible to me - a typical forced vibration
> plot that settles down to the 2.5 r/s forcing function. Is it the amplitude
> you are concerned with?
>
> --
> Ed Meyer
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://mailman.cae.wisc.edu/listinfo/help-octave
>