Uranium-238, the most prevalent isotope in uranium ore, has a half-life of about 4.5 billion years; that is, half the atoms in any sample will decay in that amount of time. See how much U238 remains for times. matlabclose all;clc;clear;% 1D radioactive decay% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 1.2 p2% Solve the Equation dN/dt = -N/tauN_uranium_initial = 1000; %initial number of uranium atomsnpoints = 1000; %Discretize time into 100 intervalsdt = 1e7; % time step in yearstau=4.4e9; % mean lifetime of 238 UN_uranium = zeros(npoints,1); % initializes N_uranium, a vector of dimension npoints X 1,to beingall zerostime = zeros(npoints,1); % this initializes the vector time to being all zerosN_uranium(1) = N_uranium_initial; % the initial condition, first entry in the vector N_uranium isN_uranium_initialtime(1) = 0; %Initialise timefor step=1:npoints-1 % loop over the timesteps and calculate the numerical solutionN_uranium(step+1) = N_uranium(step) - (N_uranium(step)/tau)*dt;time(step+1) = time(step) + dt;end% For comparison , calculate analytical solutiont=0:1e8:10e9;N_analytical=N_uranium_initial*exp(-t/tau);% Plot both numerical and analytical solution plot(time,N_uranium+5,'r',t,N_analytical,'b'); %plots the solution set(gca,'fontsize',14)xlabel('Time in years')ylabel('Number of atoms')legend('Theorical','calculate')title('Real value is shifted up 10 to compare') Note that the analytical and numerical solution are coincident in this diagram. So The Real value is shifted up 10 to compare.

Here is the code for the numerical solution of the equations of motion for a simple pendulum using the Euler method. Note the oscillations grow with time. !!
%% Euler calculation of motion of simple pendulum% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.1%clear;length= 1; %pendulum length in metresg=9.8 % acceleration due to gravitynpoints = 250; %Discretize time into 250 intervalsdt = 0.04; % time step in secondsomega = zeros(npoints,1); % initializes omega, a vector of dimension npoints X 1,to being all zerostheta = zeros(npoints,1); % initializes theta, a vector of dimension npoints X 1,to being all zerostime = zeros(npoints,1); % this initializes the vector time to being all zerostheta(1)=0.2; % you need to have some initial displacement, otherwise the pendulum willnot swingfor step = 1:npoints-1 % loop over the timestepsomega(step+1) = omega(step) - (g/length)*theta(step)*dt;theta(step+1) = theta(step)+omega(step)*dttime(step+1) = time(step) + dt;endplot(time,theta,'r' ); %plots the numerical solution in redxlabel('time (seconds) ');ylabel('theta (radians)'); ## 3.1.1 Solution using the Euler-Cromer method.

This problem with growing oscillations is addressed by performing the solution using the Euler - Cromer method. The code is below
%% Euler_cromer calculation of motion of simple pendulum% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.1%clear;length= 1; %pendulum length in metresg=9.8; % acceleration due to gravitynpoints = 250; %Discretize time into 250 intervalsdt = 0.04; % time step in secondsomega = zeros(npoints,1); % initializes omega, a vector of dimension npoints X 1,to being all zerostheta = zeros(npoints,1); % initializes theta, a vector of dimension npoints X 1,to being all zerostime = zeros(npoints,1); % this initializes the vector time to being all zerostheta(1)=0.2; % you need to have some initial displacement, otherwise the pendulum willnot swingfor step = 1:npoints-1 % loop over the timestepsomega(step+1) = omega(step) - (g/length)*theta(step)*dt;theta(step+1) = theta(step)+omega(step+1)*dt; %note that% this line is the only change between% this program and the standard Euler methodtime(step+1) = time(step) + dt;end;plot(time,theta,'r' ); %plots the numerical solution in redxlabel('time (seconds) ');ylabel('theta (radians)'); ## 3.1.2 Simple Harmonic motion example using a variety of numerical approaches

In this example I use a variety of approaches in order to solve the following, very simple, equation of motion. It is based on Equation 3.9, with k and α =1. $\frac{d^{2} y}{d t^{2}}=-y$ We take 4 approaches to solving the equation, illustrating the use of the Euler, Euler Cromer, Second order Runge-Kutta and finally the built in MATLAB® solver ODE23. The solution using the built in MATLAB® solver ODE23 is somewhat less straightforward than those using the other techniques. A discussion of the technique follows. The first step is to take the second order ODE equation and split it into 2 first order ODE equations. These are $\frac{d y}{d t}=v \quad \frac{d v}{d t}=-y$ Next you create a MATLAB® function that describes your system of differential equations. You get back a vector of times, T, and a matrix Y that has the values of each variable in your system of equations over the times in the time vector. Each column of Y is a different variable. MATLAB® has a very specific way to define a differential equation, as a function that takes one vector of variables in the differential equation, plus a time vector, as an argument and returns the derivative of that vector. The only way that MATLAB® keeps track of which variable is which inside the vector is the order you choose to use the variables in. You define your differential equations based on that ordering of variables in the vector, you define your initial conditions in the same order, and the columns of your answer are also in that order. In order to do this, you create a state vector y. Let element 1 be the vertical displacement, y1, and element 2 is the velocity,v. Next, we write down the state equations, dy1 and dy2. These are $d y 1=v$ ${d y 2=-y 1}$ Next, we create a vector dy, with 2 elements, dy1 and dy2. Finally we call the MATLAB® ODE solver ODE23. We take the output of the function called my_shm. We perform the calculation for time values range from 0 to 100. The initial velocity is 0, the initial displacement is 10. The code to do this is here
[t,y]=ode45(@my_shm,[0,100],[0,10]);

Finally, we need to plot the second column of the y matrix, containing the displacement against time. The code to do this is
plot(t,y(:,2),'r');

Here is code in all
%% Simple harmonic motion - comparison of Euler, Euler Cromer% and 2nd order Runge Kutta and built in MATLAB Runge Kutta% function ODE45to solve ODEs.% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.1% Equation is d2y/dt2 = -y% Calculate the numerical solution using Euler method in red[time,y] = SHM_Euler (10);subplot(2,2,1);plot(time,y,'r' );axis([0 100 -100 100]);xlabel('Time');ylabel('Displacement');legend ('Euler method');% Calculate the numerical solution using Euler Cromer method in blue[time,y] = SHM_Euler_Cromer (10);subplot(2,2,2);plot(time,y,'b' );axis([0 100 -20 20]);xlabel('Time');ylabel('Displacement');legend ('Euler Cromer method');% Calculate the numerical solution using Second order Runge-Kutta method in green[time,y] = SHM_Runge_Kutta (10);subplot(2,2,3);plot(time,y,'g' );axis([0 100 -20 20]);xlabel('Time');ylabel('Displacement');legend ('Runge-Kutta method');% Use the built in MATLAB ODE45 solver to solve the ODE% The function describing the SHM equations is called my_shm% The time values range from 0 to 100% The initial velocity is 0, the initial displacement is 10[t,y]=ode23(@SHM_ODE45_function,[0,100],[0,10]);% We need to plot the second column of the y matrix, containing the% displacement against time in blacksubplot(2,2,4);plot(t,y(:,2),'k');axis([0 100 -20 20]);xlabel('Time');ylabel('Displacement');legend ('ODE45 Solver');%% Simple harmonic motion - Euler method% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.1% Equation is d2y/dt2 = -yfunction [time,y] = SHM_Euler (initial_displacement);npoints = 2500; %Discretize time into 250 intervalsdt = 0.04; % time step in secondsv = zeros(npoints,1); % initializes v, a vector of dimension npoints X 1,to being all zerosy = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zerostime = zeros(npoints,1); % this initializes the vector time to being all zerosy(1)=initial_displacement; % need some initial displacement% Euler solutionfor step = 1:npoints-1 % loop over the timestepsv(step+1) = v(step) - y(step)*dt;y(step+1) = y(step)+v(step)*dt;time(step+1) = time(step) + dt;endend%% Simple harmonic motion - Euler Cromer method% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.1% Equation is d2y/dt2 = -yfunction [time,y] = SHM_Euler_Cromer (initial_displacement);npoints = 2500; %Discretize time into 250 intervalsdt = 0.04; % time step in secondsv = zeros(npoints,1); % initializes v, a vector of dimension npoints X 1,to being all zerosy = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zerostime = zeros(npoints,1); % this initializes the vector time to being all zerosy(1)=initial_displacement; % need some initial displacement% Euler Cromer solutionfor step = 1:npoints-1 % loop over the timestepsv(step+1) = v(step) - y(step)*dt;y(step+1) = y(step)+v(step+1)*dt;time(step+1) = time(step) + dt;endend%% Simple harmonic motion - Second order Runge Kutta method% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.1% Equation is d2y/dt2 = -yfunction [time,y] = SHM_Runge_Kutta(initial_displacement);% 2nd order Runge Kutta solutionnpoints = 2500; %Discretize time into 250 intervalsdt = 0.04; % time step in secondsv = zeros(npoints,1); % initializes v, a vector of dimension npoints X 1,to being all zerosy = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zerostime = zeros(npoints,1); % this initializes the vector time to being all zerosy(1)=initial_displacement; % need some initial displacementv = zeros(npoints,1); % initializes v, a vector of dimension npoints X 1,to being all zerosy = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zerosv_dash = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zerosy_dash = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zerostime = zeros(npoints,1); % this initializes the vector time to being all zerosy(1)=10; % need some initial displacementfor step = 1:npoints-1 % loop over the timesteps v_dash=v(step)-0.5*y(step)*dt; y_dash=y(step)+0.5*v(step)*dt; v(step+1) = v(step)-y_dash*dt; y(step+1) = y(step)+v_dash*dt; time(step+1) = time(step)+dt;endend%% Simple harmonic motion - Built in MATLAB ODE45 method% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.1% Equation is d2y/dt2 = -yfunction dy = SHM_ODE45_function(t,y)% y is the state vectory1 = y(1); % y1 is displacementv = y(2); % y2 is velocity% write down the state equationsdy1=v;dy2=-y1;% collect the equations into a column vector, velocity in column 1,% displacement in column 2dy = [dy1;dy2];end This solution uses q=1
%% Euler_cromer calculation of motion of simple pendulum with damping% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.2%clear;length= 1; %pendulum length in metresg=9.8; % acceleration due to gravityq=1; % damping strengthnpoints = 250; %Discretize time into 250 intervalsdt = 0.04; % time step in secondsomega = zeros(npoints,1); % initializes omega, a vector of dimension npoints X 1,to being all zerostheta = zeros(npoints,1); % initializes theta, a vector of dimension npoints X 1,to being all zerostime = zeros(npoints,1); % this initializes the vector time to being all zerostheta(1)=0.2; % you need to have some initial displacement, otherwise the pendulum willnot swingfor step = 1:npoints-1 % loop over the timestepsomega(step+1) = omega(step) - (g/length)*theta(step)*dt-q*omega(step)*dt;theta(step+1) = theta(step)+omega(step+1)*dt;% In the Euler method, , the previous value of omega% and the previous value of theta are used to calculate the new values of omega and theta.% In the Euler Cromer method, the previous value of omega% and the previous value of theta are used to calculate the the new value% of omega. However, the NEW value of omega is used to calculate the new% theta%time(step+1) = time(step) + dt;end;plot(time,theta,'r' ); %plots the numerical solution in redxlabel('time (seconds) ');ylabel('theta (radians)'); All of the next five plots were produced using the code below with slight modifications in either the input parameters or the plots.
% Euler Cromer Solution for non-linear, damped, driven pendulum% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi,% section 3.3%clear;length= 9.8; %pendulum length in metresg=9.8; % acceleration due to gravityq=0.5;F_Drive=1.2; % damping strengthOmega_D=2/3;npoints =15000; %Discretize timedt = 0.04; % time step in secondsomega = zeros(npoints,1); % initializes omega, a vector of dimension npoints X 1,to being all zerostheta = zeros(npoints,1); % initializes theta, a vector of dimension npoints X 1,to being all zerostime = zeros(npoints,1); % this initializes the vector time to being all zerostheta(1)=0.2; % you need to have some initial displacement, otherwise the pendulum willnot swingomega(1)=0;for step = 1:npoints-1;% loop over the timesteps% Note error in book in Equation for Example 3.3omega(step+1)=omega(step)+(-(g/length)*sin(theta(step))-q*omega(step)+F_Drive*sin(Omega_D*time(step)))*dt;temporary_theta_step_plus_1 = theta(step)+omega(step+1)*dt;% We need to adjust theta after each iteration so as to keep it between +/-pi% The pendulum can now swing right around the pivot, corresponding to theta>2*pi.% Theta is an angular variable so values of theta that differ by 2*pi correspond to the same position.% For plotting purposes it is nice to keep (-pi% So, if theta is <-pi, add 2*pi.If theta is > pi, subtract 2*pi% If the lines below between the ****** are commented out you get 3.6 (b)% bottom %******************************************************************************************** if (temporary_theta_step_plus_1 < -pi) temporary_theta_step_plus_1= temporary_theta_step_plus_1+2*pi; elseif (temporary_theta_step_plus_1 > pi) temporary_theta_step_plus_1= temporary_theta_step_plus_1-2*pi; end;%******************************************************************************************** % Update theta array theta(step+1)=temporary_theta_step_plus_1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In the Euler method, , the previous value of omega% and the previous value of theta are used to calculate the new values of omega and theta. % In the Euler Cromer method, the previous value of omega% and the previous value of theta are used to calculate the the new value% of omega. However, the NEW value of omega is used to calculate the new% theta%time(step+1) = time(step) + dt;end;plot (theta,omega,'r' ); %plots the numerical solutionxlabel('theta (radians)');ylabel('omega (seconds)'); Note that this code seems has problems. " % Program to perform Euler_cromer calculation of motion of physical pendulum % by Kevin Berwick, and calculate the bifurcation diagram. You need to % have the function ‘pendulum_function' available in order to run this. % based on ‘Computational Physics' book by N Giordano and H Nakanishi, % section 3.4. Omega_D=2⁄3; for F_Drive_step=1:0.1:13

F_Drive=1.35+F_Drive_step/100;
% Calculate the plot of theta as a function of time for the current drive step
% using the function :- pendulum_function
[time,theta]= pendulum_function(F_Drive, Omega_D);
%Filter the results to exclude initial transient of 300 periods, note
% that the period is 3*pi.
I=find (time< 3*pi*300);
time(I)=NaN;
theta(I)=NaN;
%Further filter the results so that only results in phase with the driving force
% F_Drive are displayed.
% Replace all those values NOT in phase with NaN
Z=find(abs(rem(time, 2*pi/Omega_D)) > 0.01);
time(Z)=NaN;
theta(Z)=NaN;
% Remove all NaN values from the array to reduce dataset size
time(isnan(time)) = [];
theta(isnan(theta)) = [];
% Now plot the results
plot(F_Drive,theta,'k');
hold on;
axis([1.35 1.5 1 3]);
xlabel('F Drive');


end % Euler_cromer calculation of motion of physical pendulum % by Kevin Berwick, % based on ‘Computational Physics' book by N Giordano and H Nakanishi, % section 3.3 function [time,theta] = pendulum_function(F_Drive,Omega_D) length= 9.8; %pendulum length in metres g=9.8; % acceleration due to gravity q=0.5; % damping strength npoints =100000; %Discretize time dt = 0.04; % time step in seconds omega = zeros(npoints,1); % initializes omega, a vector of dimension npoints X 1,to being all zeros theta = zeros(npoints,1); % initializes theta, a vector of dimension npoints X 1,to being all zeros time = zeros(npoints,1); % this initializes the vector time to being all zeros theta(1)=0.2; % you need to have some initial displacement, otherwise the pendulum will not swing omega(1)=0; for step = 1:npoints-1

% loop over the timesteps
omega(step+1)=omega(step)+(-(g/length)*sin(theta(step))-q*omega(step)+F_Drive*sin(Omega_D*time(step)))*dt;
temporary_theta_step_plus_1 = theta(step)+omega(step+1)*dt;
% Make corrections to keep theta between pi and -pi
if (temporary_theta_step_plus_1 < -pi)
temporary_theta_step_plus_1= temporary_theta_step_plus_1+2*pi;
elseif (temporary_theta_step_plus_1 > pi)
temporary_theta_step_plus_1= temporary_theta_step_plus_1-2*pi;
end
% Update theta array
theta(step+1)=temporary_theta_step_plus_1;
% Increment time
time(step+1) = time(step) + dt;


end end "

• 4.1
• 4.2
• 4.3
• 4.4
• 4.6 4.1 Kepler's Laws Johannes Kepler, working with data painstakingly collected by Tycho Brahe without the aid of a telescope, developed three laws which described the motion of the planets across the sky.
1. The Law of Orbits: All planets move in elliptical orbits, with the sun at one focus.
2. The Law of Areas: A line that connects a planet to the sun sweeps out equal areas in equal times.
3. The Law of Periods: The square of the period of any planet is proportional to the cube of the semimajor axis of its orbit. Kepler's laws were derived for orbits around the sun, but they apply to satellite orbits as well.
%% Planetary orbit using Euler Cromer methods.% by Kevin Berwick,% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 4.1%npoints=500;dt = 0.002; % time step in yearsx=1; % initialise position of planet in AUy=0;v_x=0; % initialise velocity of planet in AU/yrv_y=2*pi;% Plot the Sun at the originplot(0,0,'oy','MarkerSize',30, 'MarkerFaceColor','yellow');axis([-1 1 -1 1]);xlabel('x(AU)');ylabel('y(AU)');hold on;for step = 1:npoints-1;% loop over the timestepsradius=sqrt(x^2+y^2);% Compute new velocities in the x and y directionsv_x_new=v_x - (4*pi^2*x*dt)/(radius^3);v_y_new=v_y - (4*pi^2*y*dt)/(radius^3);% Euler Cromer Step - update positions using newly calculated velocitiesx_new=x+v_x_new*dt;y_new=y+v_y_new*dt;% Plot planet position immediatelyplot(x_new,y_new,'-k');drawnow;% Update x and y velocities with new velocitiesv_x=v_x_new;v_y=v_y_new;% Update x and y with new positionsx=x_new;y=y_new;end; (The time step should be less than 1% of the characteristic time scale of the problem.) 4.3 Precession of the perihelion of Mercury Let's do the physics here, $F{G}=\frac{G M{S} M{e}}{r^{2}}\left(1+\frac{\alpha}{r^{2}}\right)$ so $\frac{d^{2} x}{d t^{2}}=\frac{F{G, x}}{M{E}} \quad \frac{d^{2} y}{d t^{2}}=\frac{F{G, y}}{M{E}}$ Now, write each 2nd order differential equations as two, first order, differential equations. $\frac{d v{x}}{d t}=-\frac{G M{s} x}{r^{3}}\left(1+\frac{\alpha}{r^{2}}\right)$ $\frac{d x}{d t}=v{x}$ $\frac{d v{y}}{d t}=-\frac{G M{s} y}{r^{3}}\left(1+\frac{\alpha}{r^{2}}\right)$ $\frac{d y}{d t}=v{y}$ So, the difference equation set using the Euler Cromer method is $v{x, i+1}=v{x, i}-\frac{4 \pi^{2} x{i}}{r{i}^{a}} \Delta t\left(1+\frac{\alpha}{r{i}^{2}}\right)$ $x{i+1}=x{i}+v{x, i+1} \Delta t$ $v{y, i+1}=v{y, i}-\frac{4 \pi^{2} y{i}}{r{i}^{s}} \Delta t\left(1+\frac{\alpha}{r{i}^{2}}\right)$ $y{i+1}=y{i}+v{y, i+1} \Delta t$ We could go ahead and code this, but what about if we chose to attack the problem using the Runge Kutta method. The relevant equations are $y^{\prime}=y{i}+v{y, i} \frac{\Delta t}{2}$ $v{y}^{\prime}=v{y, i}-\frac{4 \pi^{2} y{i}}{r{i}^{3}} \frac{\Delta t}{2}\left(1+\frac{\alpha}{r{i}^{2}}\right)$ $y{i+1}=y{i}+v{y}^{\prime} \Delta t$ $v{y, i+1}=v{y, i}-\frac{4 \pi^{2} y{i}}{r{i}^{8}} \Delta t\left(1+\frac{\alpha}{r{i}^{2}}\right)$ and $x^{\prime}=x{i}+v{x, i} \frac{\Delta t}{2}$ $v{x}^{\prime}=v{x, i}-\frac{4 \pi^{2} x{i} \Delta t}{r{i}^{s}} \frac{\Delta t}{2}\left(1+\frac{\alpha}{r{i}^{2}}\right)$ $x{i+1}=x{i}+v{x}^{\prime} \Delta t$ $v{x, i+1}=v{x, i}-\frac{4 \pi^{2} x{i}}{r{i}^{8}} \Delta t\left(1+\frac{\alpha}{r_{i}^{2}}\right)$ 4.4 The three body problem and the effect of Jupiter on Earth A couple of points on this problem. Firstly, the direction of the various forces between the 3 bodies is elegantly captured in the pseudocode given in Ex 4.2. Do not be tempted to start taking absolute values of the subtracted positions in a naïve bid to ‘correct' the equations. Here is the code
%% 3 body simulation of Jupiter, Earth and Sun. Use Euler Cromer method% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 4.4% by Kevin Berwick%%npoints=1000000;dt = 0.0001; % time step in years.M_s=2e30; % Mass of the Sun in kgM_e=6e24; % Mass of the Earth in kgM_j=1.9e27*1e3; % Mass of Jupiter in kgx_e_initial=1; % Initial position of Earth in AUy_e_initial=0;v_e_x_initial=0; % Initial velocity of Earth in AU/yrv_e_y_initial=2*pi;x_j_initial=5.2; % Initial position of Jupiter in AU, assume at opposition initiallyy_j_initial=0;v_j_x_initial=0; % Initial velocity of Jupiter in AU/yrv_j_y_initial= 2.7549; % This is 2*pi*5.2 AU/11.85 years = 2.75 AU/year% Create arrays to store position and velocity of Earthx_e=zeros(npoints,1);y_e=zeros(npoints,1);v_e_x=zeros(npoints,1);v_e_y=zeros(npoints,1);% Create arrays to store position and velocity of Jupiterx_j=zeros(npoints,1);y_j=zeros(npoints,1);v_j_x=zeros(npoints,1);v_j_y=zeros(npoints,1);r_e=zeros(npoints,1);r_j=zeros(npoints,1);r_e_j=zeros(npoints,1);% Initialise positions and velocities of Earth and Jupiterx_e(1)=x_e_initial;y_e(1)=y_e_initial;v_e_x(1)=v_e_x_initial;v_e_y(1)=v_e_y_initial;x_j(1)=x_j_initial;y_j(1)=y_j_initial;v_j_x(1)=v_j_x_initial;v_j_y(1)=v_j_y_initial;for i = 1:npoints-1; % loop over the timesteps% Calculate distances to Earth from Sun, Jupiter from Sun and Jupiter% to Earth for current value of ir_e(i)=sqrt(x_e(i)^2+y_e(i)^2);r_j(i)=sqrt(x_j(i)^2+y_j(i)^2);r_e_j(i)=sqrt((x_e(i)-x_j(i))^2 +(y_e(i)-y_j(i))^2);% Compute x and y components for new velocity of Earthv_e_x(i+1)=v_e_x(i)-4*pi^2*x_e(i)*dt/r_e(i)^3-4*pi^2*(M_j/M_s)*(x_e(i)-x_j(i))*dt/r_e_j(i)^3;v_e_y(i+1)=v_e_y(i)-4*pi^2*y_e(i)*dt/r_e(i)^3-4*pi^2*(M_j/M_s)*(y_e(i)-y_j(i))*dt/r_e_j(i)^3;% Compute x and y components for new velocity of Jupiterv_j_x(i+1)=v_j_x(i)-4*pi^2*x_j(i)*dt/r_j(i)^3-4*pi^2*(M_j/M_s)*(x_j(i)-x_e(i))*dt/r_e_j(i)^3;v_j_y(i+1)=v_j_y(i)-4*pi^2*y_j(i)*dt/r_j(i)^3-4*pi^2*(M_j/M_s)*(y_j(i)-y_e(i))*dt/r_e_j(i)^3;%% Use Euler Cromer technique to calculate the new positions of Earth and% Jupiter. Note the use of the NEW vlaue of velocity in both equationsx_e(i+1)=x_e(i)+v_e_x(i+1)*dt;y_e(i+1)=y_e(i)+v_e_y(i+1)*dt;x_j(i+1)=x_j(i)+v_j_x(i+1)*dt;y_j(i+1)=y_j(i)+v_j_y(i+1)*dt;end;plot(x_e,y_e, 'r', x_j,y_j, 'k');axis([-7 7 -7 7]);xlabel('x(AU)');ylabel('y(AU)');title('3 body simulation - Jupiter Earth'); The 3 body simulation when Jupiter Earth is 1000 times of actual value. This 3body problem is quite interesting. 4.6 Chaotic tumbling of Hyperion The model of the moon consists of two particles, m1 and m2 joined by a massless rod. This orbits around a massive object, Saturn, at the origin. We need to extend our original planetary motion program to include the rotation of the object. First we need to recall the maths in Section 4.1, we used in order to calculate the motion of the Earth around the Sun. $\frac{d^{2} x}{d t^{2}}=\frac{F{G, x}}{M{E}} \quad \frac{d^{2} y}{d t^{2}}=\frac{F{G, y}}{M{E}}$ ${l}{F{G, x}=-\frac{G M{s} M{e}}{r^{2}} \cos \theta}$ ${F{G, x}=-\frac{G M{s} M{e} x}{r^{3}}}$ ${F{G, y}=-\frac{G M{s} M{e}}{r^{2}} \sin \theta}$ ${F{G, y}=-\frac{G M{s} M{e} y}{r^{3}}}$ Now, write each 2nd order differential equations as two, first order, differential equations ${c}{\frac{d v{x}}{d t}=-\frac{G M{s} x}{r^{3}}}$ ${\frac{d x}{d t}=v{x}}$ ${\frac{d v{y}}{d t}=-\frac{G M{s} y}{r^{3}}}$ ${\frac{d y}{d t}=v{y}}$ We need suitable units of mass. Not that the Earth's orbit is circular. For circular motion we know that the centripetal force is given by $\frac{M{E} v^{2}}{r}$ , where v is the velocity of the Earth. ${c}{\frac{M{E} v^{2}}{r}=F{G}=\frac{G M{s} M{E}}{r^{2}}}$ ${G M{s}=v^{2} r=4 \pi^{2} A U^{3} / y r^{2}}$ So, the difference equation set using the Euler Cromer method is ${l}{v{x, i+1}=v{x, i}-\frac{4 \pi^{2} x{i}}{r{i}^{3}} \Delta t}$ ${x{i+1}=x{i}+v{x, i+1} \Delta t}$ ${v{y, i+1}=v{y, i}-\frac{4 \pi^{2} y{i}}{r{i}^{3}} \Delta t}$ ${y{i+1}=y{i}+v{y, i+1} \Delta t}$ We can use this equation set to model the motion of the centre of mass of Hyperion. Now, from the analysis of the motion of Hyperion. $\omega=\frac{d \theta}{d t}$ $\frac{d \omega}{d t} \approx-\frac{3 G M{s a t}}{r{c}^{5}}\left(x{c} \sin \theta-y{c} \cos \theta\right)\left(x{c} \cos \theta+y{c} \sin \theta\right)$ So we need to add two more difference equations to our program, and noting that $G M{s a t}=4 \pi^{2}$ as noted in the book, ${c}{\omega{i+1}=\omega{i}-\frac{3\left(4 \pi^{2}\right)}{r{c}^{5}}\left(x{i} \sin \theta{i}-y{i} \cos \theta{i}\right)\left(x{i} \cos \theta{i}+y{i} \sin \theta{i}\right) \Delta t}$ ${\theta{i+1}=\theta{i}+\omega_{i+1} \Delta t}$ Here is the code for the motion of Hyperion. The initial velocity in the y direction was 1 HU/Hyperion year as explained in the book. This gave a circular orbit. Note from the results that the tumbling is not chaotic under these conditions. If we change the initial velocity in the y direction to 5 HU/Hyperion year as explained in the book. This gave an elliptical orbit. Here is the new code and below this code are the results from running this code. Note that now the motion is chaotic.
%% Simulation of chaotic tumbing of Hyperion, the moon of Saturn . Use Euler Cromer method% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 4.6% by Kevin Berwick%%npoints=100000;dt = 0.0001; % time step in yearstime=zeros(npoints,1);r_c=zeros(npoints,1);% Create arrays to store position, velocity and angle and angular velocity of% centre of massx=zeros(npoints,1);y=zeros(npoints,1);v_x=zeros(npoints,1);v_y=zeros(npoints,1);theta=zeros(npoints,1);omega=zeros(npoints,1);x(1)=1; % initialise position of centre of mass of Hyperion in HUy(1)=0;v_x(1)=0; % initialise velocity of centre of mass of Hyperionv_y(1)=5;% initialise theta and omega of Hyperionfor i= 1:npoints-1;% loop over the timestepsr_c(i)=sqrt(x(i)^2+y(i)^2);% Compute new velocities in the x and y directionsv_x(i+1)=v_x(i) - (4*pi^2*x(i)*dt)/(r_c(i)^3);v_y(i+1)=v_y(i) - (4*pi^2*y(i)*dt)/(r_c(i)^3);% Euler Cromer Step - update positions of centre of mass of Hyperion using% NEWLY calculated velocitiesx(i+1)=x(i)+v_x(i+1)*dt;y(i+1)=y(i)+v_y(i+1)*dt;% Calculate new angular velocity omega and angle theta. Note that GMsaturn=4*pi^2, see book for detailsTerm1=3*4*pi^2/(r_c(i)^5);Term2=x(i)*sin(theta(i))- y(i)*cos(theta(i));Term3=x(i)*cos(theta(i)) +y(i)*sin(theta(i));omega(i+1)=omega(i) -Term1*Term2*Term3*dt;%Theta is an angular variable so values of theta that differ by 2*pi correspond to the same position.%We need to adjust theta after each iteration so as to keep it between%+/-pi for plotting purposes. We do that heretemporary_theta_i_plus_1= theta(i)+omega(i+1)*dt;if (temporary_theta_i_plus_1 < -pi)temporary_theta_i_plus_1= temporary_theta_i_plus_1+2*pi;elseif (temporary_theta_i_plus_1 > pi)temporary_theta_i_plus_1= temporary_theta_i_plus_1-2*pi;end;% Update theta arraytheta(i+1)=temporary_theta_i_plus_1;time(i+1)=time(i)+dt;end;subplot(2,1,1);plot(time, theta,'-g');axis([0 10 -4 4]);xlabel('time(year)');ylabel('theta(radians)');title('theta versus time for Hyperion');subplot(2,1,2);plot(time, omega,'-k');axis([0 10 -20 60]);xlabel('time(year)');ylabel('omega(radians/yr)');title('omega versus time for Hyperion'); 5. Potentials and Fields 5.1 Solution of Laplace's equation using the Jacobi relaxation method. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Jacobi method to solve Laplace equation% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 5.1% by Kevin Berwick%% Load array into V[V] =Initialise_V_Jacobi_metal_box;% run update routine and estimate convergence%Initialise loop counterloops=1;[V_new, delta_V_new]=Update_V_Jacobi_Metal_box(V);% While we have not met the convergence criterion and the number of loops is <10 so that we give the relaxation% algorithm time to convergewhile (delta_V_new > 49e-5 & loops < 10);loops=loops+1;[V_new, delta_V_new]=Update_V_Jacobi_Metal_box(V_new);% draw the surface using the mesh functionmesh (V_new);title('Potential Surface');drawnow;% insert a pause here so we see the evolution of the potential% surfacepause(1);end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Jacobi method to solve Laplace equation% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 5.1% by Kevin Berwick%% This function creates the intial voltage array Vfunction[V] =Initialise_V_Jacobi_metal_box;% clear variablesclear;V = [-1 -0.67 -0.33 0 0.33 0.67 1;-1 0 0 0 0 0 1;-1 0 0 0 0 0 1;-1 0 0 0 0 0 1;-1 0 0 0 0 0 1;-1 0 0 0 0 0 1;-1 -0.67 -0.33 0 0.33 0.67 1];end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function [ V_new, delta_V_new] = Update_V_Jacobi_Metal_box(V);% This function takes a matrix V and applies Eq 5.10 to it. Only the values% inside the boundaries are changed. It returns the % processed matrix to% the calling function, together with the value of delta_V, the total% accumulated amount by which the elements% of the matrix have changedrow_size = size(V,1);column_size=size(V,2);% preallocate memory for speedV_new=V;delta_V_new=0;% Move along the matrix, element by element computing Eq 5.10, ignoring% boundariesfor j =2:column_size-1;for i=2:row_size -1;V_new(i,j) = (V(i-1,j)+V(i+1,j)+V(i,j-1)+V(i,j+1))/4;% Calculate delta_V_new value, the cumulative change in V during this% update call, to test for convergencedelta_V_new=delta_V_new+abs(V_new(i,j)-V(i,j));end;end;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5.1.1 Solution of Laplace's equation for a hollow metallic prism with a solid, metallic inner conductor. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Jacobi method to solve Laplace equation% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 5.1% by Kevin Berwick%% Load array into V[V] =Initialise_prism;% run update routine and estimate convergence[V_new, delta_V_new]=Update_prism(V);%Initialise loop counterloops=0;% While we have not met the convergence criterion and the number of loops% is <10 so that we give the relaxation% algorithm time to converge% while (delta_V_new > & loops < 30);while (delta_V_new>4e-5 | loops < 20);loops=loops+1;[V_new, delta_V_new]=Update_prism(V_new);% draw the surface using the mesh function% mesh (V_new,'FaceColor','interp','EdgeColor','none','FaceLighting','phong');mesh (V_new,'FaceColor','interp');title('Potential Surface');axis([0 20 0 20 0 1]);drawnow;% insert a pause here so we see the evolution of the potential% surfacepause(0.5);end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Jacobi method to solve Laplace equation% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 5.1% by Kevin Berwick%% This function creates the intial voltage array Vfunction[V] =Initialise_prism;% clear variablesclear;V = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 00 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 00 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 00 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 00 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 00 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function [ V_new, delta_V_new] = Update_prism(V);% This function takes a matrix V and applies Eq 5.10 to it. Only the values% inside the boundaries are changed. It returns the processed matrix to% the calling function, together with the value of delta_V, the total% accumulated amount by which the elements % of the matrix have changedrow_size = size(V,1);column_size=size(V,2);% preallocate memeory for speedV_new=V;delta_V_new=0;% Move along the matrix, element by element computing Eq 5.10, ignoring% boundariesfor j =2:column_size-1;for i=2:row_size -1;% Do not update V in metal barif V(i,j) ~=1;% If the value of V is not =1, calculate V_new and% delta_V_new to test for convergenceV_new(i,j) = (V(i-1,j)+V(i+1,j)+V(i,j-1)+V(i,j+1))/4;delta_V_new=delta_V_new+abs(V_new(i,j)-V(i,j))else% otherwise, leave value unchangedV_new(i,j)=V(i,j);end;end;end; end 5.1.2 Solution of Laplace's equation for a finite sized capacitor %% Jacobi method to solve Laplace equation% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 5.1% by Kevin Berwick%% Load array into V[V] =capacitor_initialise;% run update routine and estimate convergence[V_new, delta_V_new]=capacitor_update(V);%Initialise loop counterloops=1;% While we have not met the convergence criterion and the number of loops% is <20 so that we give the relaxation% algorithm time to convergewhile (delta_V_new>400e-5 | loops < 20);loops=loops+1;[V_new, delta_V_new]=capacitor_update(V_new);% draw the surface using the mesh functionmesh (V_new,'Facecolor','interp');title('Potential Surface');axis([0 20 0 20 -1 1]);drawnow;% insert a pause here so we see the evolution of the potential% surfacepause(0.5);end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Jacobi method to solve Laplace equation% based on 'Computational Physics' book by N Giordano and H Nakanishi% Section 5.1% by Kevin Berwick%% This function creates the intial voltage array Vfunction[V] =capacitor_initialise;% clear variablesclear;V = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 00 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function [ V_new, delta_V_new] = capacitor_update(V);% This function takes a matrix V and applies Eq 5.10 to it. Only the values% inside the boundaries are changed. It returns the % processed matrix to% the calling function, together with the value of delta_V, the total% accumulated amount by which the elements % of the matrix have changed row_size = size(V,1);column_size=size(V,2);% preallocate memory for speedV_new=zeros(row_size, column_size);delta_V_new=0;% Move along the matrix, element by element computing Eq 5.10, ignoring% boundariesfor j =2:column_size-1;for i=2:row_size -1;% Do not update V on the platesif V(i,j)~=1 & V(i,j) ~=-1;% If the value of V is not =1 or -1, calculate V_new and% cumulative delta_V_new to test for convergenceV_new(i,j) = (V(i-1,j)+V(i+1,j)+V(i,j-1)+V(i,j+1))/4;delta_V_new=delta_V_new+abs(V_new(i,j)-V(i,j))else% otherwise, leave value unchangedV_new(i,j)=V(i,j);end;end;end;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% code of this course