Experiences with Scilab – 3
Let us try differential equations. We want to solve the following problem.
\(y“ = -g + k |y’|^a, \quad y(0)=h.\)
This simulates a falling body with friction. We want to know, when the body hits the ground. After reformulating the problem as a system of differential equations, the following code in Euler Math Toolbox does the job.
>function f(x,y) := [y;-g+k*abs(y)^a] >g=9.81; k=0.01; a=1.67; h=100; >x=0:0.01:8; plot2d(x,ode("f",x,[h;0]),xl="time",yl="height"): >function g(x) := ode("f",[0,x],[h;0])[1,2]; >solve("g",5) 4.78765346143
EMT uses the LSODA method. This is an adaptive method. So we need no intermediate points between time 0 and time x in the call to ode(). Note that Since I have no doubts about the quality of the routines in Scilab this should be easy in Scilab too. However, as a newbie it took some time to figure out the mistakes I made. Here is the result.
function yd = ydot(x,y) yd = [y(2);-g+k*abs(y(2))^a] endfunction g=9.81; k=0.01; a=1.67; h=100; t=0:0.01:8; y1=ode([h;0],0,t,ydot) plot(t,y1(1,:)) function res = yh(x) res=ode([h;0],0,x,ydot) res=res(1) endfunction y1=fsolve(5,yh)
I still do not know how to print something from a script. So I store y1 to a variable and print it afterwards.
-->format(16) -->y1 y1 = 4.787653545896
The solutions disagree somewhat. I checked that with other methods. EMT is closer to the solution. However, there might be options in Scilab I do not know.
The mistakes I made at first try are numerous. Of course, I tend to write y instead of y(2) for vectors. But it also took a long time to discover the error when I defined the function as
function y = ydot(x,y) ydot = [y(2);-g+k*abs(y(2))^a] endfunction
I also discovered the somewhat disappointing fact that res=ode(…)(1) is not possible in the function yh().
I observed some differences between the systems.
- Functions can see global variables in Scilab. In EMT, only one-line functions or functions with the „useglobal“ command can see global variables.
- If the function ydot() needed extra parameters Scilab uses a list to pack the function name and the extra parameters. This is a good idea and a clean solution. In EMT, you have to use semicolon parameters.
- The function ode() can call a C routine ydot() in Scilab. The routine can be in C form. In EMT, the C routine must handle the stack of EMT. I will explore C routines in Scilab on some other day.
- Functions and variables with the same name cannot coexist. Probably, that is a good idea for a system which uses round brackets to get elements of vectors. In EMT, this is not so necessary.