Reputation: 29
i've been trying to solve a second order differential equation using this 1, however I couldnt get it right and found nothing helpful online but I believe I've made progress.
I used dsolve;
syms x(t) v(t) fi(t)
[x(t), v(t)] = dsolve(diff(x) == v, diff(v) == fi/m, x(0) == [-L, -L], v(0) == [5, 10] )
which gives me;
x(t) =
int(fi(x)/5, x, 0, t, 'IgnoreSpecialCases', true, 'IgnoreAnalyticConstraints', true) - 5
v(t) =
C2 + t*(int(fi(x)/5, x, 0, t, 'IgnoreSpecialCases', true, 'IgnoreAnalyticConstraints', true) - 5) + int(-(x*fi(x))/5, x, 0, t, 'IgnoreSpecialCases', true, 'IgnoreAnalyticConstraints', true)
Now I need help interperating the result and I wonder if I can use this result to get something out of ode45? Also I want to plot the solution as a reference trajectory to a simulation of 500 particles moving through a forcefield.
using ODE45:
function dxdt = solution(t,y0)
frprintf('Second stop')
.....
dxdt = [x, v]
end
call from main file:
t:dt:t_f
y0 = [x0,v0]
fprintf('first stop')
[x, v] = ode45(@solution, y0, t)
I've set it up so that if the code would run smoothly it would print "first stop, second stop, third stop" and "fourth stop", it only prints out first stop and thats where I get the error.
Upvotes: 0
Views: 1317
Reputation: 25972
Assuming that you somewhere defined
function F = fi(t)
F = ...
end
you define the ODE function as (using m
as global variable)
function doty = odefunc(t,y)
doty = [ y(2); fi(t)/m ]
end
and call
t = t0:dt:tf
y0 = [ x0, v0 ]
t,y = ode45(odefunc, t, y0)
plot( t, y(:,1) )
Generally, the vector y
would contain a point in the phase space (position and velocity/impulse) for all particles or objects in the system. In odefunc
you would then compute the forces for this particular phase space point at the time t
and compose from that the derivative vector to the phase space point.
In a 3D simulation you, for instance, could arrange that y(6*(k-1)+1:6*(k-1)+3)
are the position of particle k
and y(6*(k-1)+4:6*(k-1)+6)
the velocity vector. Or you could separate position and velocity with y(3*(k-1)+1:3*(k-1)+3)
for position and y(3*(N+k-1)+1:3*(N+k-1)+3)
for the velocity for the k
-th of N
particles.
Upvotes: 1