# 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[2];-g+k*abs(y[2])^a]
>g=9.81; k=0.01; a=1.67; h=100;
>x=0:0.01:8; plot2d(x,ode("f",x,[h;0])[1],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[2] 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.

Diese Website verwendet Akismet, um Spam zu reduzieren. Erfahre mehr darüber, wie deine Kommentardaten verarbeitet werden.