lsoda-- caution...
Walter Lehmann
I wrote a scilab script to simulate a swinging system. If I reduce the
moment of inertia (in the script the variable "J" is used) I receive the
message:
lsoda-- caution... t (=r1) and h (=r2) are
such that t + h = t at next step
(h = pas). integration continues
where r1 is : 0.1946504419246D-02 and r2 :
0.1055861088740D-18
I guess, that this means that I should lower the step size of the solver.
Could anybody help me?
Thanks in advance,
Walter
The script follows:
clear
function f = deq (t, y, Ma, fqz, lc, c, J)
f(1) = y(2)
f(2) = (1 / J) * ((Ma * sin (2 * %pi * fqz * t)) - ((lc^2) * c * tan (y
(1))))
endfunction
// Amplitude des Antriebsmoments [N m]
Ma = 10;
// Antriebsfrequenz [Hz]
fqz = 50;
// Kraftangriffspunkt der Feder [m]
lc = 8E-3;
// Federkonstante [N / m]
c = 2.56E3;
// Trägheitsmoment [kg m^2]
J = 2.3697E-6;
// Simulationsdauer [s]
Te = 0.1;
// Anzahl Datenpunkte
dp = 1000;
t = linspace (0, Te, dp);
y0 = [0; 0];
y = ode (y0, 0, t, deq);
xgrid (1)
plot2d (t, y(1,:))
RRogers
I haven't run the code but adjusting rtol, atol is probably the way to
go; alternately "rfk" and "fix" procedures. %eps (D-16) is the
resolution in Scilab and you have to keep the relative step size above
that.
http://turing.une.edu.au/~amth247/Lectures_2003/Lecture_06/lecture/
has an easy discussion of computational limits; but I doubt if you are
truly being bothered by them.
If the problem is stiff then you might use the "stiff" option.
Rainer von Seggern
> has an easy discussion of computational limits; but I doubt if you are
> truly being bothered by them.
> If the problem is stiff then you might use the "stiff" option.
Hello
I think you should use a list in the call of ode. The following script
runs under Scilab 5.1.1:
function f = deq(t,y,Ma,fqz,lc,c,J)
f(1)=y(2)
f(2)=Ma*sin(2*%pi*fqz*t)/J -lc^2*c*tan(y(1))
endfunction
Ma=10; fqz=50; lc=8E-3; c=2.56E;
J=2.3697E-6; Te=0.1; dp=1000;
t=linspace(0,Te,dp);
y0=[0;0];
y=ode(y0,0,t,list(deq,Ma,fqz,lc,c,J));
xgrid(1)
plot2d(t,y(1,:))
Good luck, Rainer
Walter Lehmann
Hi Rainer,
I modified the code like you suggested above. It doesn' make any
difference.
Cheers Stefan
Walter Lehmann
Hi,
sorry for my late response. Could you give me a hint how I could adjust
these parameteres?
Thank you,
Walter
Rainer von Seggern
Hi Stefan,
I don't understand that. Here is the result of my script for y and
dp=20:
0. 0.
28.086424 14541.744
148.43226 26681.791
253.54064 10135.018
268.90625 727.81525
314.3311 18828.276
444.53425 25246.032
530.67626 6085.6225
539.31707 2832.3744
604.81899 22530.085
738.43101 22530.097
803.93303 2832.3933
812.57394 6085.6416
898.71604 25246.051
1028.9193 18828.304
1074.3443 727.84473
1089.7101 10135.049
1194.8187 26681.828
1315.1647 14541.778
1343.2513 0.0343614
Alles Gute für 2010,
Rainer
Walter Lehmann
Hello Rainer,
thank you, your right. It also works for me. Could you explain to me why
your version works, mine not.
Thanks and also a happy new year,
Walter
Rainer von Seggern
Hello Walter,
I think the reason might be, that one of the parameters Ma, fqz,...
which are defined in the calling script might be overwritten in the
Scilab-program ode.
Using the list construct this is not possible. I think I have not yet
understand all the mysteries of using functions in Scilab.
Rainer
Rainer von Seggern
Hello Walter,
I think now I have found the true reason of your problem.
First of all I was using a wrong function:
function f = deq(t,y,Ma,fqz,lc,c,J)
f(1)=y(2)
f(2)=Ma*sin(2*%pi*fqz*t)/J - lc^2*c*tan(y(1))
endfunction
while your function is
function f = deq(t,y,Ma,fqz,lc,c,J)
f(1)=y(2)
f(2)=(Ma*sin(2*%pi*fqz*t) - lc^2*c*tan(y(1)))/J
endfunction
I think in your case the argument y(1) of the tangens comes near %pi/
2, and this singularity is the problem !
Try the following for t1=2e-3, and you will see what I mean:
t1=1.9e-3;y=ode(y0,0,t1,list(deq,Ma,fqz,lc,c,J))'
You have to change the model to avoid the singularity.
Good luck,
Rainer