How to solve system of 2nd order differential equations using ode45

20 views (last 30 days)
siyeong Jang on 16 Nov 2019
Answered: Stephan on 16 Nov 2019
I have three 2nd order differential equations with my initial conditions and I'm trying to use the ode45 function in matlab to solve this.
I wish get (x,y) position in x-y plane
but i can`t simultaneously get x-position and y-position respect t.
%planar circular restricted 3 body problem
% %Initial Conditions
% x(0) = 35.0; %AU
% y(0) = 0; %AU
%Constants
G = 1;
m1 = 0.52; %solar mass
m2 = 0.33; %solar mass
r1= -21.3; %AU
r2= 33.7; %AU
ome=sqrt((m1+m2)/(-r1+r2)^3);
% x=x(1) x'=x(2) y=x(3) y'=x(4)
% x''=x*ome^2+2*ome*y'-m1*(x-r1)/(((x-r1)^2+y^2)^1.5)-m2*(x-r2)/(((x-r2)^2+y^2)^1.5)
% y''=y*ome^2-2*ome*x'-m1*y/(((x-r1)^2+y^2)^1.5)-m2*y/(((x-r2)^2+y^2)^1.5)
% x=x1 , x1'=x2 . y=x3 , x3'=x4
[t,x]=ode45(@crtbp,[0 35],[35;0;0;0.03]);
function dxdt=crtbp(t,x)
m1=0.52;
m2=0.33;
r1=-21.3;
r2=33.7;
ome=sqrt((m1+m2)/(-r1+r2)^3);
dxdt=[x(2); x(1)*ome^2+2*ome*x(4)-m1*(x(1)-r1)/(((x(1)-r1)^2+x(3)^2)^1.5)-m2*(x(1)-r2)/(((x(1)-r2)^2+x(3)^2)^1.5); x(4); x(3)*ome^2-2*ome*x(2)-m1*x(3)/(((x(1)-r1)^2+x(3)^2)^1.5)-m2*x(3)/(((x(1)-r2)^2+x(3)^2)^1.5)];
end

1 Comment

darova on 16 Nov 2019
Looks correct. What is the problem? You want you plot x vs y?
plot(x(:,1),x(:,3))