Solving a System of ODEs

My Input:

syms x(t) y(t) z(t)
%Units of kg and seconds
q = 1.60217662*(10^-19);
B = 2;
m_e = 1.60217662*(10*-31);
a = (q*B)/m_e;

%DiffiQ to solve
ode1 = diff(x,2) == a*diff(y);
ode2 = diff(y,2) == -a*diff(x);
ode3 = diff(z,2) == 0;
odes = [ode1; ode2; ode3];

%Solutions
S = dsolve(odes);
xSol(t) = S.x;
ySol(t) = S.y;
zSol(t) = S.z;

[xSol(t), ySol(t), zSol(t)] = dsolve(odes);

%Initial Conditions
condx1 = x(0) == 1;
condy1 = y(0) == 1;
condz1 = z(0) == 0;

%{
condvx1 = diff(x) == 1;
condvy1 = diff(y) == 2;
condvz1 = diff(z) == 2;
%}

%Plot Trajectory

t = 0:pi/50:8*pi;
plot3(ode1,ode2,ode3,'r','LineWidth',3)

My Output: Nothing. The program runs but nothing happens. I have to force close MatLab so that the script stops. I left it running for an hour to see but I am starting to think that there is another issue. Any advice is appreciated!

PS - Using Matlab R2018a and yes I purposefully commented out a section of initial conditions.

Answer 1

[Star Strider]

Try this:

...

%Initial Conditions
condx1 = x(0) == 1;
condy1 = y(0) == 1;
condz1 = z(0) == 0;

%Solutions
S = dsolve(odes, condx1, condy1, condz1);
xSol(t) = S.x;
ySol(t) = S.y;
zSol(t) = S.z;

Sx = matlabFunction(xSol);
Sy = matlabFunction(ySol);
Sz = matlabFunction(zSol);

Then evaluate the anonymous functions in terms of the variables and constants. Note that you will have to either evaluate ‘Sz’ first, then pass the result as ‘z’ or pass it as an evaluated function ‘Sz(t,Cz)’ as the ‘z’ argument.

[Tom Keaton]

Sorry, but can you please explain the "matlabFunction" that is being used? Is this the equivalent to an "anonymous function" or is this another arbitrary command? Also to clarify, when you say pass the result, do you mean do this?:

Sx = matlabFunction(xSol);
Sy = matlabFunction(ySol);
Sz = matlabFunction(zSol);

%Pass results
x = Sx;
y = Sy;
z = Sz;

[Star Strider]

The matlabFunction (link) function creates anonymous functions (or function files) from symbolic functions. Here, the calls create:

Sx =

    function_handle with value:

      @(t,C14,z)z+cos(t.*6.451612903225807e-22)-sin(t.*6.451612903225807e-22)-z.*cos(t.*6.451612903225807e-22)+C14.*sin(t.*6.451612903225807e-22)

Sy =

    function_handle with value:

      @(t,C14,z)C14+cos(t.*6.451612903225807e-22)+sin(t.*6.451612903225807e-22)-z.*sin(t.*6.451612903225807e-22)-C14.*cos(t.*6.451612903225807e-22)

Sz =

    function_handle with value:

      @(t,C13)C13.*t

So you would evaluate and plot them as:

t = 0:pi/50:8*pi;
C13 = ...;
C14 = ...;
xv = Sx(t,C14,Sz(t,C13));
yv = Sx(t,C14,Sz(t,C13));
zv = Sz(t,C13);

plot3(xv, yv, zv, 'r','LineWidth',3)
grid on

Another option would be to evaluate ‘Sz’ first, and pass ‘zv’ as the ‘z’ argument.

[Tom Keaton]

(UPDATE) I see and thank you for all this assistance! I read through the matlabFunction text also. This is very interesting. I also tried doing a run with a stop at every line and found that the script has an error in line 9 (Aka the first coupled ODE in this system). So there is something wrong with the fact that MatLab seems unable to solve even the first ODE alone. Do you know why this may be? I posted the error message below:

Error using mupadmex
Internal error with symbolic engine. Quit and restart MATLAB.

Error in sym>tomupad (line 1219)
    S = mupadmex(numeric2cellstr(x));

Error in sym (line 211)
                S.s = tomupad(x);

Error in syms (line 201)
toDefine = sym(zeros(1, 0));

Error in parttraject (line 1)
syms x(t) y(t) z(t)

[Tom Keaton]

PS - I already quit and restarted 3 different times

[Star Strider]

My pleasure.

The error you stated and the error message you posted (thank you) appear to be completely different problems. The error appears to be in Line 1 of your code, where you first call the Symbolic Math Toolbox. If that is actually the problem, see: "syms", "sym", and "mupad" functions cause MATLAB to freeze (link). This has the appropriate solution.

Meanwhile, please post the corrected differential equations so I can see if the corrected set change the results in my code.

[Tom Keaton]

Alright, I looked at the link and just to clarify, the fix is to download the latest update? If so, I already have the latest version I think.

[v d] = version

v =

      '9.4.0.813654 (R2018a)'

d =

      'February 23, 2018'

Also I know for certain these coupled ODEs are correct based on my analytical calculations and double checking them from other resources.

[Star Strider]

Install the update anyway. It won’t change anything that doesn’t need to be changed.

I can’t help you with ODEs I don’t have.

[Walter Roberson]

You need the patches on top of that release.

[Tom Keaton]

@Star Strider Sorry for the late response. When you say "corrected ODEs" do you mean the ODEs in plain text form? If so, here they are:

[x" = a*y'], [y" = -a*x'], [z" = 0] Note: All derivatives are wrt time so for example: [x" = d^2x/dt^2] and [x' = dx/dt]

Where: [a = (q*B)/m_e] (A single, simple constant to simplify the term out front), [q = 1.60217662*(10^-19)] (Charge of electron), [B = 2] (Constant B-Field in the z direction), [m_e = 1.60217662*(10*-31)] (Mass of electron)

[Tom Keaton]

@Walter Roberson I tried looking for an update, but there was nothing on the site that prompted me to do so. Should I run the installer again? Or is there a way to do it through the program while I have it open? Like a menu option?

[Walter Roberson]

https://www.mathworks.com/downloads/web_downloads/download_update?release=R2018a&s_tid=ebrg_R2018a_1_1755918

[Star Strider]

@Tom Keaton — You mentioned that you had changed one of your differential equations. Apparently you haven’t, since everything is the same.

If you want to use the Symbolic Math Toolbox after the recent Windows 10 Update, you need to install the MATLAB update I linked to. It also updates several other Toolboxes, if you have them. Note Walter’s Comment.

Also, one option is to integrate your functions numerically, although the output is disappointing. (To do this, you need to add Y to your syms declaration.)

That aside:

syms x(t) y(t) z(t) Y
...
[VF,Subs] = odeToVectorField(odes);
odefcn = matlabFunction(VF, 'Vars',{t, Y});

t = linspace(0, 8*pi, 50);
ic = [1; 0; 1; 0; 0; 0];
[t,y] = ode15s(odefcn, t, ic);

figure
plot3(y(:,1), y(:,3), y(:,5), '-p')
grid on
xlabel('X')
ylabel('Y')
zlabel('Z')

Multiplying the constant by 1E+15 does not change the results, so the magnitude is not the problem.

[Tom Keaton]

Let me try all these the next time I get back to this simulation portion. I will respond again with results (If any) in a few days.

[Star Strider]

Noted.

[Tom Keaton]

Hey @Star Strider. I am back and will be engaged in this thread again. I talked with Matlab staff and they sent me a specialized link to get the correct update to fix the bug. The bug is now fixed, so now it is just about getting the code I have to work.

[Star Strider]

Noted.

The update in "syms", "sym", and "mupad" functions cause MATLAB to freeze worked for me.

[Tom Keaton]

I was able to get this to work now actually. I found out that Matlab's ODEs Toolbox just doesn't support systems of higher order differntial equations. It was only "recently" too that this language is able to solve higher order differential equations in the first place. So I was just forced to create 6, first order differential equations and the system was able to solve them. Here is the code:

syms x(t) y(t) z(t) vx(t) vy(t) vz(t)
%Units of kg and seconds
q = 1.60217662*(10^-19);
B = 2;
m_e = 1.60217662*(10*-31);
a = (q*B)/m_e;

%DiffiQ to solve
ode1 = diff(x) == vx;
ode2 = diff(y) == vy;
ode3 = diff(z) == vz;
ode4 = diff(vx) == a*vy;
ode5 = diff(vy) == -a*vx;
ode6 = diff(vz) == 0;
odes = [ode1; ode2; ode3; ode4; ode5; ode6];

%Initial Conditions
condx1 = x(0) == 1;
condy1 = y(0) == 1;
condz1 = z(0) == 0;
condvx1 = vx(0) == (1*10^6);
condvy1 = vy(0) == (2*10^6);
condvz1 = vz(0) == (2*10^6);
conds = [condx1; condy1; condz1; condvx1; condvy1; condvz1];

%Solutions
S = dsolve(odes,conds);
xSol(t) = S.x
ySol(t) = S.y
zSol(t) = S.z

%Plot Trajectory
%t = 0:pi/50:16*pi*m*(1/(q*B));
t = linspace(0,16*pi*m_e*(1/(q*B)),5000);

figure(1)
plot3(xSol(t),ySol(t),zSol(t),'r','LineWidth',3)

[Star Strider]

I can’t find any reference to ‘ODEs Toolbox’ in the documentation. However the ODE solvers in core MATLAB have no problems with your system:

syms x(t) y(t) z(t) vx(t) vy(t) vz(t) Y 
%Units of kg and seconds
q = 1.60217662*(10^-19);
B = 2;
m_e = 1.60217662*(10*-31);
a = (q*B)/m_e;

%DiffiQ to solve
ode1 = diff(x) == vx;
ode2 = diff(y) == vy;
ode3 = diff(z) == vz;
ode4 = diff(vx) == a*vy;
ode5 = diff(vy) == -a*vx;
ode6 = diff(vz) == 0;
odes = [ode1; ode2; ode3; ode4; ode5; ode6];

[ODEsVF, Subs] = odeToVectorField(odes);
odesfcn = matlabFunction(ODEsVF, 'Vars',{t,Y});
icv = [1; 1; 0; 1E+6; 2E+6; 2E+6];
t = linspace(0,16*pi*m_e*(1/(q*B)),5000);

[T,Y] = ode15s(odesfcn, t, icv);

figure
plot3(Y(:,2), Y(:,1), Y(:,3), '-r', 'LineWidth',3)
grid on

This uses odeToVectorField to create a vector field from your ‘odes’ array, and matlabFunction to convert it to an anonymous function that the ODE solvers can use. I use ode15s because the constants in your system vary by several orders-of-magnitude, and such systems are usually ‘stiff’. Note that the plot arguments appear to be out of order. See the ‘Subs’ variable to understand the reason.

[Tom Keaton]

I see now why you were using the MatlabFunction now. I am still trying to wrap my head around how it works which is probably why I didn't see how it would solve this issue before. I understand now though why it was used but I still need to read up more about the function and what it does exactly. So is it the case then that one may use this function to solve any higher order set of coupled and uncoupled ODEs in Matlab?

[Star Strider]

‘So is it the case then that one may use this function to solve any higher order set of coupled and uncoupled ODEs in Matlab?’

In general, yes. However the required calls are to odeToVectorField first, and then to matlabFunction.

There are other functions, such as odeFunction (link) and related functions that are useful for differential-algebraic equations (that do not apply to your system here).

With your system, pay particular attention to the ‘Subs’ variable in my code. Those values correspond to the ‘Y’ values in the odeToVectorField output, and to the order matlabFunction specifies and the integrated output columns that the ODE solvers return. This is the reason the plot3 arguments in my code appear out-of-order. They correspond to your plot3 call, and are in the order the functions return.

[Tom Keaton]

I see. The past few days I have been messing around with these and trying to implement them in other ways. I will keep messing around with them, especially since the equations I have right now are only simplified versions of what I really am trying to do and I can only apply this separation trick to these. I will close this thread and accept the answer since the original question has been answered. I will be opening up more threads in the near future as I continue developing this simulation. Thank you again for all the help and patience!

[Star Strider]

As always, my pleasure!