Hello Friends,

I am having a problem in solving the following 2nd order differential equation. Please help me out as to which method would be best to tackle it.

The Problem statement is to find the trajectory of moon, in a 3 body problem.
The code is as follows:
---------------------------------------------------------------------------------------
% ------PartI: The simulation of sun earth moon system-------
clear;
clc;

Me=5.972E24;                 % Mass of earth
Ms= 1.988E30;                % Mass of sun
R= 1.495E11;                 % Orbital radius of earth
r= 3.85E8;                   % Orbital radius of moon
G= 6.67E-11;                 % Gravitational Constant
T=3.155E7;                   % orbital time of earth
w= sqrt(G*(Me+Ms))/R;        % Angular velocity
de= Ms*R/(Ms+Me);            % Distance between earth and centre of mass
ds= Me*R/(Ms+Me);            % Distance between sun and centre of mass
% vv=1.018;                  % average velocity of moon

% Propagation of time
t=linspace(0,80,100);

% Equation for coordinates of sun and earth
xs=ds*cos(w.*t);
ys=ds*sin(w.*t);
xe=-de*cos(w.*t);
ye=-de*sin(w.*t);

% Plotting of trajectory of sun and earth
fig = figure;
ax = axes;
plot(xs,ys,'yo');
hold on;
plot(xe,ye,'go');
grid on;
set(ax,'color','black');

% -----Part II: Simulation of trajectory of moon around earth and sun-----

x(size(t))=0;
y(size(t))=0;
a= sqrt(((xs - x).^2) + ((ys - y).^2));
b= sqrt(((xe - x).^2) + ((ye - y).^2));
p= ((G*Ms.*(xs-x))./power(a,3) + (G*Me.*(xe-x))./power(b,3));
q= ((G*Ms.*(ys-y))./power(a,3) + (G*Me.*(ye-y))./power(b,3));
X=dsolve('D2x-p=0','x(0)= -(de+r)','Dx(0)=0');
Y=dsolve('D2y-q=0','y(0)=0','Dy(0)=0');
x=subs(X);
y=subs(Y);
plot(x,y,'c.');
-------------------------------------------------------------------------------------------------------

Here x, y, xs, xe, ys, ye are the functions of t.

that gives us the p and q as function of t as well. But if you see the resulting equations X and Y, you will find that they are being treated as a constant.
I fail to understand what am I doing wrong here.

Later, I also thought about using ode45 to solve the equations, by first converting them to 1st order using the following help page:
http://www.mathworks.com/support/tech-notes/1500/1510.html

But the complexity of the equation itself is not allowing me to break it into pure form, having separate x(t) and y(t) (see the denominator of q and w, which has component of x and y in there).

Unless I can somehow separate the function x and y, the ode45 method won't give me much of a result...

Please help me out on how to properly use the dsolve and/or ode45 for these equation. If there are any doubts in what the problem, do ask.

Also give me any pointers if you can, to make it a better code.

Thank You
Nishant 

"Nishant Gupta" <im_nishantg@yahoo.co.in> wrote in message 
news:ifhq71$85h$1@fred.mathworks.com...
> Hello Friends,
>
> I am having a problem in solving the following 2nd order differential 
> equation. Please help me out as to which method would be best to tackle 
> it.

*snip*

> X=dsolve('D2x-p=0','x(0)= -(de+r)','Dx(0)=0');
> Y=dsolve('D2y-q=0','y(0)=0','Dy(0)=0');

You don't appear to have any symbolic variables in your expression, so I 
would use the NUMERIC ODE solvers (like ODE45) instead of the SYMBOLIC ODE 
solver DSOLVE.

If you must use DSOLVE for this problem, substitute your function for p into 
the expression BEFORE calling DSOLVE.

fcn = subs('D2x-p=0', 'p', p);
X=dsolve(fcn,'x(0)= -(de+r)','Dx(0)=0');

*snip*

> Later, I also thought about using ode45 to solve the equations, by first 
> converting them to 1st order using the following help page:
> http://www.mathworks.com/support/tech-notes/1500/1510.html
>
> But the complexity of the equation itself is not allowing me to break it 
> into pure form, having separate x(t) and y(t) (see the denominator of q 
> and w, which has component of x and y in there).

Let z = [x; dxdt; y; dydt].  You can now write:

dzdt = [z(2); yourExpressionForD2x; z(4); yourExpressionForD2y]

which is just:

dzdt = [z(2); p; z(4); q]

Given t and z as input your ODE function should be able to calculate p and q 
just fine.  Use z(1) and z(3) in place of x and y respectively.

-- 
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlab.wikia.com/wiki/FAQ
To contact Technical Support use the Contact Us link on 
http://www.mathworks.com 

"Steven_Lord" <slord@mathworks.com> 

Thanks a lot for the help. I used ode45 in the way you suggested.
the following are the two files:

moony.m
-----------------------------------------------------------------------------------------------------
function dz = moony(t,z)
% function to be integrated

% The constants required in the calculation

Me=5.972E24;                 % Mass of earth
Ms= 1.988E30;                % Mass of sun
R= 1.495E11;                 % Orbital radius of earth
% r= 3.85E8;                 % Orbital radius of moon
G= 6.67E-11;                 % Gravitational Constant
% TT=3.155E7;                % orbital time of earth
w= sqrt(G*(Me+Ms)/R^3);      % Angular velocity
de= Ms*R/(Ms+Me);            % Distance between earth and centre of mass
ds= Me*R/(Ms+Me);            % Distance between sun and centre of mass
% vv=1.018E3;                % average velocity of moon

dz = zeros(4,1);

% Writing Equation for gravitational acceleration on moon

a= sqrt(((ds*cos(w*t) - z(1))^2) + ((ds*sin(w*t) - z(3))^2));
b= sqrt(((de*cos(w*t) + z(1))^2) + ((de*sin(w*t) + z(3))^2));
p = ((G*Ms*(ds*cos(w*t)-z(1)))/power(a,3) - (G*Me*(de*cos(w*t)+z(1)))/power(b,3));
q = ((G*Ms*(ds*sin(w*t)-z(3)))/power(a,3) - (G*Me*(de*sin(w*t)+z(3)))/power(b,3));

% Giving the values to function dz

dz(1) = z(2);
dz(2) = p;
dz(3) = z(4);
dz(4) = q;

---------------------------------------------------------------------------------------------------------------
moony2.m
---------------------------------------------------------------------------------------------------------------
% The constants required in the calculation

Me=5.972E24;                 % Mass of earth
Ms= 1.988E30;                % Mass of sun
R= 1.495E11;                 % Orbital radius of earth
r= 3.85E8;                   % Orbital radius of moon
G= 6.67E-11;                 % Gravitational Constant
TT=3.155E7;                  % orbital time of earth
w= sqrt(G*(Me+Ms)/R^3);      % Angular velocity
de= Ms*R/(Ms+Me);            % Distance between earth and centre of mass
ds= Me*R/(Ms+Me);            % Distance between sun and centre of mass
vv=1.018E3;                  % average velocity of moon

% Defining tspan, initial condition.
tspan = linspace(0,TT,1000);

z0 = [-(de+r); 0; 0; vv];
[t,z] = ode45(@moony, tspan, z0);


% Equation for coordinates of sun and earth
xsun= ds.*cos(w.*tspan);
ysun= ds.*sin(w.*tspan);
xea= -de.*cos(w.*tspan);
yea= -de.*sin(w.*tspan);

% Plotting of trajectory of sun and earth
fig = figure;
ax = axes;
plot(xsun,ysun,'yo');
hold on;
plot(xea,yea,'go');
grid on;
set(ax,'color','black');
plot(z(:,1),z(:,3),'wo-');
-----------------------------------------------------------------------------------------------------------
most of the problems are solved, but the moon is still not behaving the way it should...
I checked my equations and initial conditions, but couldn't find anything wrong yet.
Following is the crude explanation of the theory:

The sun and earth are a 2 body system. The centre of mass of sun and moon is placed at the origin. 'ds' and 'de' are the distance from the center of mass (origin) of sun and earth resp.

de +ds =R

according to center of mass formula:
de= Ms*R/(Ms+Me)
ds= Me*R/(Ms+Me)

we know that

Ms*ds*w^2 = G*Ms*Me/(R^2)   -- centripetal and centrifugal forces

from here, we get 'w' - angular velocity.

Now, gravitational force at any point (x,y) in 2 dimensional space (at any time t), due to sun and moon are given by 'p' and 'q'. as the two bodies (sun and moon) are revolving about COM, with a constant angular velocity 'w'....

http://en.wikipedia.org/wiki/Newton%27s_law_of_universal_gravitation

The initial condition:

Sun is on the positive x-axis (at distance ds). Earth is on negative x axis (at distance de). and moon is on left hand side of earth, on x axis (at distance 'r') with velocity in the upward direction....