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!