Matlab Solvers: A Case Study in Mechanics each other. For example, think of the
ID: 2084789 • Letter: M
Question
Matlab Solvers: A Case Study in Mechanics
each other. For example, think of the earth and the moon, where the moon moves
around the earth at distance 1. (Of course, here both the masses and the distance are
normalized.) A third object, which is relatively much smaller and does not affect the
motion of the first two, is also orbiting in space. Think of this as a space-ship, or a
meteorite.
To simplify the analysis, we assume that all trajectories lie in the same plane. We choose
the barycenter of mass as the origin of our coordinate system. Moreover, we adopt a
rotating frame of coordinates. In this moving frame, the earth and the moon are always
located at the points
and their barycenter is at (0,0). We denote by (y1(t), y2(t)) the position of the space ship at time t, w.r.t. this rotating frame of coordinates.
The equations of motion are:
Notice that, in each of above ODEs, the first two terms account for the Coriolis force, due
to the fact that the rotating frame is not an inertial frame. The last two terms account
for the gravitational pull of the earth and the moon, respectively. The numbers dE and dM
measure the distance of the space ship from the earth and from the moon, respectively. Consider the initial
conditions:
For this particular initial data, the solution turns out to be periodic. Indeed, at the time t = 17.065211656 the rotating object goes back to it's initial configuration, and the entire motion is repeated.
You need to do the following:
Explanation / Answer
follow the code step by step which i was written for the ODE equations
***********
*******************************
%ODE equations in matlab
function f=fun1(t,y)
f=-t*y/sqrt(2-y^2);
[tv1 f1]=ode23('fun1',[0 5],1);
[tv2 f2]=ode45('fun1',[0 5],1);
plot(tv1,f1,'-.',tv2,f2,'--')
title('y''=-ty/sqrt(2-y^2), y(0)=1, t in [0, 5]')
grid
axis([0 5 0 1])
%intial values
y1'=2y1+y2+5y3+e-2t
y2'=-3y1-2y2-8y3+2e-2t-cos(3t)
y3'=3y1+3y2+2y3+cos(3t)
y1(0)=1, y2(0)=-1, y3(0)=0
t in [0,pi/2].
%F(t,Y) for any given t, y1, y2, and y3 and name it funsys.m
function Fv=funsys(t,Y);
Fv(1,1)=2*Y(1)+Y(2)+5*Y(3)+exp(-2*t);
Fv(2,1)=-3*Y(1)-2*Y(2)-8*Y(3)+2*exp(-2*t)-cos(3*t);
Fv(3,1)=3*Y(1)+3*Y(2)+2*Y(3)+cos(3*t);
***************************
%Eulers method
h=0.1; % step's size
N=2400; % number of steps
y(1)=1;
for n=1:N
y(n+1)= y(n)+h*(-6*y(n));
x(n+1)=n*h;
end
plot(x,y)
***************************
%4th order Runge kutta method
clc
clf
clear all
% Solve dy/dt = f(t,y). In general it can be a function of both
% define f(t,y)
fcnstr='y-t^2+1' ;
f=inline(fcnstr,'t','y') ;
% t0, initial time
t0=0 ;
% y0, corresponding value of y at t0
y0=0.5;
% tf, upper boundary of time t (end time).
tf=2;
% n, number of steps to take
n=5;
%**********************************************************************
% Displays title
disp(sprintf(' The 4th Order Runge-Kutta Method of Solving Ordinary Differential Equations'))
h=(tf-t0)/n ;
disp(sprintf(' h = ( tf - t0 ) / n '))
disp(sprintf(' = ( %g - %g ) / %g ',tf,t0,n))
disp(sprintf(' = %g',h))
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
disp(sprintf(' Step %g',i))
disp(sprintf('---------------------------'))
% Adding Step Size
ta(i+1)=ta(i)+h ;
% Calculating k1, k2, k3, and k4
k1 = f(ta(i),ya(i)) ;
k2 = f(ta(i)+0.5*h,ya(i)+0.5*k1*h) ;
k3 = f(ta(i)+0.5*h,ya(i)+0.5*k2*h) ;
k4 = f(ta(i)+h,ya(i)+k3*h) ;
% Using 4th Order Runge-Kutta formula
ya(i+1)=ya(i)+1/6*(k1+2*k2+2*k3+k4)*h ;
disp('1) Find k1 and k2 using the previous step information.')
disp(sprintf(' k1 = f( x%g , y%g )',i-1,i-1))
disp(sprintf(' = f( %g , %g )',ta(i),ya(i)))
disp(sprintf(' = %g ',k1))
disp(sprintf(' k2 = f( x%g + 0.5 * h , y%g + 0.5 * k1 * h )',i-1,i-1))
disp(sprintf(' = f( %g + 0.5 * %g , %g + 0.5 * %g * %g)',ta(i),h,ya(i),k1,h))
disp(sprintf(' = f( %g , %g )',ta(i)+0.5*h,ya(i)+0.5*k1*h))
disp(sprintf(' = %g ',k2))
disp(sprintf(' k3 = f( x%g + 0.5 * h , y%g + 0.5 * k2 * h )',i-1,i-1))
disp(sprintf(' = f( %g + 0.5 * %g , %g + 0.5 * %g * %g)',ta(i),h,ya(i),k2,h))
disp(sprintf(' = f( %g , %g )',ta(i)+0.5*h,ya(i)+0.5*k2*h))
disp(sprintf(' = %g ',k3))
disp(sprintf(' k4 = f( x%g + h, y%g + k3*h)',i-1,i-1))
disp(sprintf(' = f( %g , %g )',ta(i)+h,ya(i)+k3*h))
disp(sprintf(' = %g ',k1))
disp(sprintf('2) Apply the Runge-Kutta 4th Order method to estimate y%g',i))
disp(sprintf(' y%g = y%g + 1/6*( k1 + 2*k2 + 2*k3 +k4) * h',i,i-1))
disp(sprintf(' = %g + %g * %g * %g',ya(i),1/6,(k1+2*k2+2*k3+k4),h))
disp(sprintf(' = %g ',ya(i+1)))
disp(sprintf(' at t%g = %g',i,ta(i+1)))
end
disp(sprintf(' **********************Results*******************'))
% The following finds what is called the 'Exact' solution
tspan = [t0 tf];
[t,y]=ode45(f,tspan,y0);
[yfi dummy]=size(y);
yf=y(yfi);
% Plotting the Exact and Approximate solution of the ODE.
hold on
xlabel('x');ylabel('y');
title('Exact and Approximate Solution of the ODE by the 4th Order Runge-Kutta Method');
plot(t,y,'--','LineWidth',2,'Color',[0 0 1]);
plot(ta,ya,'-','LineWidth',2,'Color',[0 1 0]);
legend('Exact','Approximation');
disp('The figure window that now appears shows the approximate solution ')
disp('piecewise continuous straight lines The blue line represent the exact')
disp('solution. In this case ''exact'' refers to the solution obtained by the')
disp(sprintf('Matlab function ode45. '))
disp('Note the approximate value obtained as well as the true value and ')
disp('relative true error at our desired point x = xf.')
disp(sprintf(' Approximate = %g',ya(n+1)))
disp(sprintf(' Exact = %g',yf))
disp(sprintf(' True Error = Exact - Approximate'))
disp(sprintf(' = %g - %g',yf,ya(n+1)))
disp(sprintf(' = %g',yf-ya(n+1)))
disp(sprintf(' Absolute Relative True Error Percentage'))
disp(sprintf(' = | ( Exact - Approximate ) / Exact | * 100'))
disp(sprintf(' = | %g / %g | * 100',yf-ya(n+1),yf))
disp(sprintf(' = %g',abs( (yf-ya(n+1))/yf )*100))
disp(sprintf(' The approximation can be more accurate if we made our step'))
disp(sprintf('size smaller (that is, increasing the number of steps). '))
********************************