Re: lsode precision

[Juan Pablo Carbajal]

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

jwe answered in the bug report and he saw the error on the original
code. You have to plot
y(:,1) not y(:,2)
the second signal is the velocity.

I am sorry for noticing that right away.

In any case, always check vs. analytic solutions.

And as H E Dudeney said
"We have to use our best judgment as to whether a puzzle contains a
catch or not; but we should never hastily assume it. To quibble over
the conditions is the last resort of the defeated would-be solver."

Cheers