solving a system of 4 second order differential equations

2 次查看（过去 30 天）
Rebecca Roper
Rebecca Roper 2020-3-23
评论： Rena Berman 2020-5-14
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 个评论
显示 3更早的评论
Stephen23
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
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

更多回答（1 个）

Steven Lord
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.

类别

在 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!

Translated by