이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
Solving a System of ODEs
조회 수: 1 (최근 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 (한국어)