You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Solving a System of ODEs
5 views (last 30 days)
Show older comments
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.
Accepted Answer
Star Strider
on 20 Jun 2018
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.
22 Comments
Tom Keaton
on 20 Jun 2018
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
on 21 Jun 2018
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
on 21 Jun 2018
Edited: Tom Keaton
on 21 Jun 2018
(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
on 21 Jun 2018
PS - I already quit and restarted 3 different times
Star Strider
on 21 Jun 2018
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
on 21 Jun 2018
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
on 21 Jun 2018
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
on 21 Jun 2018
You need the patches on top of that release.
Tom Keaton
on 23 Jun 2018
Edited: Tom Keaton
on 23 Jun 2018
@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
on 23 Jun 2018
@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
on 23 Jun 2018
Star Strider
on 23 Jun 2018
@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
on 26 Jun 2018
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.
Tom Keaton
on 12 Jul 2018
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.
Tom Keaton
on 12 Jul 2018
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
on 12 Jul 2018
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
on 16 Jul 2018
Edited: Tom Keaton
on 16 Jul 2018
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
on 16 Jul 2018
‘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
on 19 Jul 2018
Edited: Tom Keaton
on 19 Jul 2018
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
on 19 Jul 2018
As always, my pleasure!
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)