ODE45 in 3 dimensions with 6 initial conditions

How can you get a solution for an ODE45 problem that involves x,y,z position and velocity initial condtions and has second order equations of motion. I can't find any help with anything that is 3-dimensional and has 6 componenets for the initial conditions.

Answer 1

[James Tursa]

I find it easier to keep the position elements together, and the velocity elements together, in the state vector. So I would define the state vector y this way:

y(1) = x position

y(2) = y position

y(3) = z position

y(4) = x velocity

y(5) = y velocity

y(6) = z velocity

Then you write a derivative function based on that state vector. It would look something like this:

function ydot = myderivative(t,y,param1,param2) % you fill in the actual parameters you need
pos = y(1:3);
vel = y(4:6);
acc = _____; % you fill this in for acceleration. An expression involving pos, vel, param1, param2
ydot = [vel;acc];
end

The acc expression in the above code would be based on your 2nd order derivative expression for the actual problem you are working. I just used generic names param1 and param2 above. But you might have actual parameters named mu, g, Cd, k, etc. Then your main code would look something like this:

param1 = whatever; % or whatever the actual parameters are for your particular problem
param2 = whatever;
pos0 = whatever; % a 3x1 vector
vel0 = whatever; % a 3x1 vector
t0 = whatever;
tend = whatever;
[t,y] = ode45(@(t,y)myderivative(t,y,param1,param2),[t0 tend],[pos0;vel0]);
% code here to examine y or plot y etc.

Give this a shot and if you still have problems don't hesitate to follow up here with more questions.

[Tommy McDonald]

I got it to work! Thank you so much for the help, I was completley lost. Just didn't know how to bring in the third component. Thank you!

[Elif]

Hello, I used your solution for my code of calculating a trajectory of an electron in varying electric and magnetic fields. For some reason, when I set the starting position and velocities to zero, the code does not calculate the electric or magnetic fields when they are position dependent.

Below is my derivative function

function Va = deriv(t,rV,m,q)
% derivative function to calculate trajectory in ode45 solver

r = rV(1:3); %position vector
V = rV(4:6); %velocity vector

E = field_E(r,t); %electric field
B = field_B(r,t); %magnetic field
crs = cross(V,B); %VxB

a = [(q/m)*(E(1)+crs(1));(q/m)*(E(2)+crs(2));(q/m)*(E(3)+crs(3))];

Va = [V;a];

end

[James Tursa]

You would need to show us your equations for field_E( ) and field_B( ) to verify, but I suspect they may be inverse square law based on position. So if the position is 0 then the calculation divides by 0. Check to see if that is the case.