help-octave
[Top][All Lists]
Advanced

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

Re: lsode precision


From: Juan Pablo Carbajal
Subject: Re: lsode precision
Date: Tue, 22 Jan 2013 13:10:04 +0100

On Tue, Jan 22, 2013 at 1:09 PM, Juan Pablo Carbajal
<address@hidden> wrote:
> 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
>>

I am sorry, I posted my answer at the top of the message!


reply via email to

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