Solving a System of ODEs
3 次查看(过去 30 天)
显示 更早的评论
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.
采纳的回答
Star Strider
2018-6-20
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 个评论
Tom Keaton
2018-6-20
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
2018-6-21
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
2018-6-21
编辑:Tom Keaton
2018-6-21
(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)
Star Strider
2018-6-21
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
2018-6-21
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
2018-6-21
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.
Tom Keaton
2018-6-23
编辑:Tom Keaton
2018-6-23
@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
2018-6-23
@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
2018-6-23
Star Strider
2018-6-23
@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
2018-6-26
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
2018-7-12
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
2018-7-12
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
2018-7-12
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
2018-7-16
编辑:Tom Keaton
2018-7-16
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
2018-7-16
‘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
2018-7-19
编辑:Tom Keaton
2018-7-19
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!
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)