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

dgl

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.

 

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht.

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