How to pull specific data from solution set?

Hello, I am using this ODE solver and I am trying to figure out how I can reference the velocity components of the matrix of solutions. I think I would know how to pull the position components given the assigned variables below because the matrix "W" contains them, but how may I be able to define a velocity matrix of solutions? For context, I am solving the Lorentz force equations. The x, y and z and vx, vy and vz components are predefined and can be completely random. The "vparts" and "sparts" are the velocity and position components of arrays.

Main Script:

[X,Y,Z] = meshgrid(-0.5:.01:0.5,-0.5:.01:0.5,-0.5:.01:0.5);
[Bx, By, Bz] = B_test();
Bfieldx = arrayfun(Bx,X,Y,Z);
Bfieldy = arrayfun(By,X,Y,Z);
Bfieldz = arrayfun(Bz,X,Y,Z);

icv = [x; y; z; vx; vy; vz];

%Time Span (sec)
tspan = [0 tstep];
[T,W] = ode15s(@bdipuniodefun, tspan, icv);
[rownum,colnum] = size(W);

plot3(W(:,1), W(:,2), W(:,3), '-r', 'LineWidth',2,'color',[randi(0:1) randi(0:1) randi(0:1)])
xlabel 'x';
ylabel 'y';
zlabel 'z';
grid on
      
%Redfine the velocity and position components to reference on next if-loop run
vparts(1) = vx;
vparts(2) = vy;
vparts(3) = vz;
sparts(1) = W(rownum,1);
sparts(2) = W(rownum,2);
sparts(3) = W(rownum,3);

B_Test:

function [Bx, By, Bz] = B_test()

syms x y z
mu_0_red = 10E-7;
m = [0,0,1.28];
r = [x, y, z];
B = mu_0_red *(((dot(m,r)*r*3)/norm(r)^5) - m/norm(r)^3);
Bx = matlabFunction(B(1));
By = matlabFunction(B(2));
Bz = matlabFunction(B(3));

bdipolefun:

function bdip = bdipuniodefun(t,s)

%Using SI units
q = 1.60217662E-19;
m_e = 9.11E-31;
[Bx, By, Bz] = B_test();
bdip = [s(4); s(5); s(6); (q/m_e)*(s(5)*Bz(s(1),s(2),s(3)) - s(6)*By(s(1),s(2),s(3))); (q/m_e)*(s(6)*Bx(s(1),s(2),s(3)) - s(4)*Bz(s(1),s(2),s(3))); (q/m_e)*(s(4)*By(s(1),s(2),s(3)) - s(5)*Bx(s(1),s(2),s(3)))];

Answer 1

[Bruno Luong]

What prevent you to use bdipuniodefun() to compute the velocity dW/dt from T and W ?

[Tom Keaton]

Well the function file only computes the positions and that is what "W" represents. I am not sure just how to go about calculating DW/dt. Modify the function file "bdipolefun"? Or the original definition of [T,W]?

[Bruno Luong]

What do you mean? bdipuniodefun(t,W) must return dW/dt, that's how ode solvers requires to get in order to integrate the solution W(t).

Or I'm missing something????

[Tom Keaton]

The odes15 (link) solver returns time and position matrices, not derivatives of the functions. The coupled differential equations involved in the calculation involve both second and first derivatives so it's not as trivial as just integrating them.

[Bruno Luong]

Yeah the doc in this link one can read:

"The function dydt = odefun(t,y), for a scalar t and a column vector y, must return a column vector dydt of data type single or double that corresponds to f(t,y)."

So your function bdipuniodefun(t,W) returns dW/dt.

Simply use it to get the derivative.

You don't know what you are coding?

The term "integrating" is use for ODE not only for integration.

[Bruno Luong]

Obviously the position is column 1,2,3 of W

P = W(:,1:3);

So

dWdt = bdipuniodefun(t,W);
% velocity
dPdt = dWdt(:,1:3); % which is == W(:,4:6) according to his code of bdipuniodefun

So simply

     velocity = W(:,4:6)

You seem to just get confused big time and can't see something too obvious right there.

[Tom Keaton]

Yes, I am not too certain what this ode solver is doing which is why I am confused. So the next set of columns 4-6 of the "W" matrix are the velocity components then?

[Bruno Luong]

Evidently yes

[Tom Keaton]

I just manually checked each component and this is true. I wasn't sure how this output of solutions looked so thank you for clarifying that.