[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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
-------------------------------------------------------------