Oskar Paulsson
Oskar Paulsson

Reputation: 29

Interpereting dsolve output to use with ODE45 [MATLAB]

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

Answers (1)

Lutz Lehmann
Lutz Lehmann

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

Related Questions