help-octave
[Top][All Lists]

## Re: How can I solve equation of motion?

 From: Przemek Klosowski Subject: Re: How can I solve equation of motion? Date: Mon, 2 Dec 2002 10:00:48 -0600

```
Now,I research to solve the equation of motion by using octave.
The following is the example of equations.

(dx/dt)''=F - (dx/dt)'

#"F" is the externalforce which is measured by experiment.This value is
refered to data of experiment every time.

I came up with the program to solve this equation as follows.
But,this program does not seem to work correctly.

F=rand(5,1); # Now, I give randam "F"
function dx = ex(x,t)
global n;
dx(1) = F(n+1,1)-x(1);
dx(2) = x(1);
endfunction

global n;
v=0; # set initial value
X=0; # set initial value
for n=0:4
x0 = [v;X]   # initial value
t = linspace(n,n+1,2)
y = lsode("ex",x0,t)
v=(2,1)              # new initial value
X=(2,2)              # new initial value
end

Traditionally, the Newton equation of motion is written as F = m x''
Your equation has inconsistent dimensions (the units on the left are
m/s^3, while the right side subtracts m/s^2 from kg*m/s^2 (assuming
that you define your force in

You would solve F=mx'' by converting it into a vector first order
equation: X=[x,x']; then X'=[x',x''] = [x', F/m] = [0 1; 0 0] * X + [0 F/m].

X0=[0 0];  # X0 = [x,v];
T=linspace(0,10,100);
function xp = ex(x,t)
xp=[0 1 ; 0 0]*x + [0 3]'; # F/m = 3, for instance
endfunction
r=lsode("ex",X0,T);

gives what I think is the expected result (x = 3/2 t^2, x' = 3t). Note
that you don't need to split up the interval and solve separately---you
can just modify F in the function F, depending on the interval.

-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

```