solving a system of 4 second order differential equations
2 ビュー (過去 30 日間)
古いコメントを表示
i have the following first order equations
z1=x, z2=dx/dt, z3=y and z4=dy/dt
which relate to the following second order equations
x'=z2,x''=4*z(4)+z(1)-((M*(z(1)+E))/re^3)-((E*(z(1)-M))/rm^3), y'=z4 and y''=2*z(2)+z(3)-((M*z(3))/re^3)-((E*z(3))/rm^3)]
i have created the following function file to slove the system using ode45. however i get the error index exceeds number of array elements.
function zprime=secode(t,z)
%secode: Computes the derivatives involved in solving the first order
%system
zprime=[z(2);2*z(4)+z(1)-((M*(z(1)+E))/re^3)-((E*(z(1)-M))/rm^3);z(4);-2*z(2)+z(3)-((M*z(3))/re^3)-((E*z(3))/rm^3)];
this is my solution file code
xspan=input('please enter your time range:');
initp1=input('please enter initial condition for x(0)=?:');
initp2=input('please enter initial condition for y(0)=?:');
initv1=input('please enter the initial condition for dx/dt(0)=?:');
initv2=input('please enter the initial condition for dy/dt(0)=?:');
cond1=[initp1 initp2];
cond2=[initv1 initv2];
[x y]=ode45(@secode,xspan,cond1,cond2);
re=sqrt((z(1)+E)^2+z(3)^2);
rm=sqrt((z(1)-M)^2+z(3)^2);
M=1-E;
E=0.0123;
any help would be greatly appriciated.
5 件のコメント
Stephen23
2020 年 3 月 24 日
Original question from Google Cache (and formatted properly):
"solving a system of 4 second order differential equations"
i have the following first order equations
z1=x, z2=dx/dt, z3=y and z4=dy/dt
which relate to the following second order equations
x'=z2,x''=4*z(4)+z(1)-((M*(z(1)+E))/re^3)-((E*(z(1)-M))/rm^3), y'=z4 and y''=2*z(2)+z(3)-((M*z(3))/re^3)-((E*z(3))/rm^3)]
i have created the following function file to slove the system using ode45. however i get the error index exceeds number of array elements.
function zprime=secode(t,z)
%secode: Computes the derivatives involved in solving the first order
%system
zprime=[z(2);2*z(4)+z(1)-((M*(z(1)+E))/re^3)-((E*(z(1)-M))/rm^3);z(4);-2*z(2)+z(3)-((M*z(3))/re^3)-((E*z(3))/rm^3)];
this is my solution file code
xspan=input('please enter your time range:');
initp1=input('please enter initial condition for x(0)=?:');
initp2=input('please enter initial condition for y(0)=?:');
initv1=input('please enter the initial condition for dx/dt(0)=?:');
initv2=input('please enter the initial condition for dy/dt(0)=?:');
cond1=[initp1 initp2];
cond2=[initv1 initv2];
[x y]=ode45(@secode,xspan,cond1,cond2);
re=sqrt((z(1)+E)^2+z(3)^2);
rm=sqrt((z(1)-M)^2+z(3)^2);
M=1-E;
E=0.0123;
any help would be greatly appriciated.
採用された回答
darova
2020 年 3 月 23 日
I believe that re and rm should be functions too
re = @(z) sqrt((z(1)+E)^2+z(3)^2);
rm = @(z) sqrt((z(1)-M)^2+z(3)^2);
zprime = @(t,z) [z(2)
2*z(4)+z(1)-((M*(z(1)+E))/re(z)^3)-((E*(z(1)-M))/rm(z)^3)
z(4)
-2*z(2)+z(3)-((M*z(3))/re(z)^3)-((E*z(3))/rm(z)^3)];
tspan = [0 5e-3];
z0 = [0 1 0 1];
[t, z] = ode45(zprime,tspan,z0);
Pay attention
0 件のコメント
その他の回答 (1 件)
Steven Lord
2020 年 3 月 23 日
The third input argument to ode45 is the vector of initial conditions and the fourth is the options. You're trying to pass vectors of initial conditions in as the third and fourth inputs, and by looking at the third input ode45 thinks you only have two ODEs not four. That means it's going to call your ODE function with a two-element vector z. Combine the initial condition vectors into one vector and pass the combined vector as the third input.
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!