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. matlab
close 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/tau
N_uranium_initial = 1000; %initial number of uranium atoms
npoints = 1000; %Discretize time into 100 intervals
dt = 1e7; % time step in years
tau=4.4e9; % mean lifetime of 238 U
N_uranium = zeros(npoints,1); % initializes N_uranium, a vector of dimension npoints X 1,to being
all zeros
time = zeros(npoints,1); % this initializes the vector time to being all zeros
N_uranium(1) = N_uranium_initial; % the initial condition, first entry in the vector N_uranium is
N_uranium_initial
time(1) = 0; %Initialise time
for step=1:npoints-1 % loop over the timesteps and calculate the numerical solution
N_uranium(step+1) = N_uranium(step) - (N_uranium(step)/tau)*dt;
time(step+1) = time(step) + dt;
end
% For comparison , calculate analytical solution
t=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')

udeclay 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 metres
g=9.8 % acceleration due to gravity
npoints = 250; %Discretize time into 250 intervals
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
for step = 1:npoints-1 % loop over the timesteps
omega(step+1) = omega(step) - (g/length)*theta(step)*dt;
theta(step+1) = theta(step)+omega(step)*dt
time(step+1) = time(step) + dt;
end

plot(time,theta,'r' ); %plots the numerical solution in red
xlabel('time (seconds) ');
ylabel('theta (radians)');

chap31

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 metres
g=9.8; % acceleration due to gravity
npoints = 250; %Discretize time into 250 intervals
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
for step = 1:npoints-1 % loop over the timesteps
omega(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 method
time(step+1) = time(step) + dt;
end;
plot(time,theta,'r' ); %plots the numerical solution in red
xlabel('time (seconds) ');
ylabel('theta (radians)');
chap311

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 black
subplot(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 = -y
function [time,y] = SHM_Euler (initial_displacement);
npoints = 2500; %Discretize time into 250 intervals
dt = 0.04; % time step in seconds
v = zeros(npoints,1); % initializes v, a vector of dimension npoints X 1,to being all zeros
y = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zeros
time = zeros(npoints,1); % this initializes the vector time to being all zeros
y(1)=initial_displacement; % need some initial displacement
% Euler solution
for step = 1:npoints-1 % loop over the timesteps
v(step+1) = v(step) - y(step)*dt;
y(step+1) = y(step)+v(step)*dt;
time(step+1) = time(step) + dt;
end
end


%
% 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 = -y
function [time,y] = SHM_Euler_Cromer (initial_displacement);
npoints = 2500; %Discretize time into 250 intervals
dt = 0.04; % time step in seconds
v = zeros(npoints,1); % initializes v, a vector of dimension npoints X 1,to being all zeros
y = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zeros
time = zeros(npoints,1); % this initializes the vector time to being all zeros
y(1)=initial_displacement; % need some initial displacement
% Euler Cromer solution
for step = 1:npoints-1 % loop over the timesteps
v(step+1) = v(step) - y(step)*dt;
y(step+1) = y(step)+v(step+1)*dt;
time(step+1) = time(step) + dt;
end
end
%
% 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 = -y
function [time,y] = SHM_Runge_Kutta(initial_displacement);
% 2nd order Runge Kutta solution
npoints = 2500; %Discretize time into 250 intervals
dt = 0.04; % time step in seconds
v = zeros(npoints,1); % initializes v, a vector of dimension npoints X 1,to being all zeros
y = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zeros
time = zeros(npoints,1); % this initializes the vector time to being all zeros
y(1)=initial_displacement; % need some initial displacement
v = zeros(npoints,1); % initializes v, a vector of dimension npoints X 1,to being all zeros
y = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zeros
v_dash = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zeros
y_dash = zeros(npoints,1); % initializes y, a vector of dimension npoints X 1,to being all zeros
time = zeros(npoints,1); % this initializes the vector time to being all zeros
y(1)=10; % need some initial displacement
for 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;
end
end

%
% 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 = -y
function dy = SHM_ODE45_function(t,y)
% y is the state vector
y1 = y(1); % y1 is displacement
v = y(2); % y2 is velocity
% write down the state equations
dy1=v;
dy2=-y1;
% collect the equations into a column vector, velocity in column 1,
% displacement in column 2
dy = [dy1;dy2];
end

chap312

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 metres
g=9.8; % acceleration due to gravity
q=1; % damping strength
npoints = 250; %Discretize time into 250 intervals
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
for step = 1:npoints-1 % loop over the timesteps
omega(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 red
xlabel('time (seconds) ');
ylabel('theta (radians)');
chap32

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 metres
g=9.8; % acceleration due to gravity
q=0.5;
F_Drive=1.2; % damping strength
Omega_D=2/3;
npoints =15000; %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
% Note error in book in Equation for Example 3.3
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;
% 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 solution
xlabel('theta (radians)');
ylabel('omega (seconds)');
chap33

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');
ylabel('theta (radians)');

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 "`

  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 years
    x=1; % initialise position of planet in AU
    y=0;
    v_x=0; % initialise velocity of planet in AU/yr
    v_y=2*pi;
    % Plot the Sun at the origin
    plot(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 timesteps
    radius=sqrt(x^2+y^2);
    % Compute new velocities in the x and y directions
    v_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 velocities
    x_new=x+v_x_new*dt;
    y_new=y+v_y_new*dt;
    % Plot planet position immediately
    plot(x_new,y_new,'-k');
    drawnow;
    % Update x and y velocities with new velocities
    v_x=v_x_new;
    v_y=v_y_new;
    % Update x and y with new positions
    x=x_new;
    y=y_new;
    end;
    chap41 (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) $ chap42 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 kg
    M_e=6e24; % Mass of the Earth in kg
    M_j=1.9e27*1e3; % Mass of Jupiter in kg
    x_e_initial=1; % Initial position of Earth in AU
    y_e_initial=0;
    v_e_x_initial=0; % Initial velocity of Earth in AU/yr
    v_e_y_initial=2*pi;
    x_j_initial=5.2; % Initial position of Jupiter in AU, assume at opposition initially
    y_j_initial=0;
    v_j_x_initial=0; % Initial velocity of Jupiter in AU/yr
    v_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 Earth
    x_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 Jupiter
    x_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 Jupiter
    x_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 i

    r_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 Earth
    v_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 Jupiter
    v_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 equations
    x_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');

    chap44 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. chap45 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 years
    time=zeros(npoints,1);
    r_c=zeros(npoints,1);
    % Create arrays to store position, velocity and angle and angular velocity of
    % centre of mass
    x=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 HU
    y(1)=0;
    v_x(1)=0; % initialise velocity of centre of mass of Hyperion
    v_y(1)=5;
    % initialise theta and omega of Hyperion
    for i= 1:npoints-1;

    % loop over the timesteps
    r_c(i)=sqrt(x(i)^2+y(i)^2);
    % Compute new velocities in the x and y directions
    v_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 velocities
    x(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 details
    Term1=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 here
    temporary_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 array
    theta(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');

    chap46 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 counter
    loops=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 converge
    while (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 function
    mesh (V_new);
    title('Potential Surface');
    drawnow;
    % insert a pause here so we see the evolution of the potential
    % surface
    pause(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 V
    function[V] =Initialise_V_Jacobi_metal_box;
    % clear variables
    clear;
    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 changed
    row_size = size(V,1);
    column_size=size(V,2);
    % preallocate memory for speed
    V_new=V;
    delta_V_new=0;
    % Move along the matrix, element by element computing Eq 5.10, ignoring
    % boundaries
    for 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 convergence

    delta_V_new=delta_V_new+abs(V_new(i,j)-V(i,j));
    end;
    end;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    chap51 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 counter
    loops=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
    % surface
    pause(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 V
    function[V] =Initialise_prism;

    % clear variables
    clear;
    V = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 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 changed
    row_size = size(V,1);
    column_size=size(V,2);
    % preallocate memeory for speed
    V_new=V;
    delta_V_new=0;
    % Move along the matrix, element by element computing Eq 5.10, ignoring
    % boundaries
    for j =2:column_size-1;
    for i=2:row_size -1;

    % Do not update V in metal bar
    if V(i,j) ~=1;
    % If the value of V is not =1, calculate V_new and
    % delta_V_new to test for convergence
    V_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 unchanged
    V_new(i,j)=V(i,j);
    end;
    end;
    end;
    end
    chap511 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 counter
    loops=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 converge
    while (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 function
    mesh (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
    % surface
    pause(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 V
    function[V] =capacitor_initialise;
    % clear variables
    clear;
    V = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 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 speed
    V_new=zeros(row_size, column_size);
    delta_V_new=0;
    % Move along the matrix, element by element computing Eq 5.10, ignoring
    % boundaries
    for j =2:column_size-1;
    for i=2:row_size -1;

    % Do not update V on the plates
    if 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 convergence
    V_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 unchanged
    V_new(i,j)=V(i,j);
    end;
    end;
    end;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    chap511 code of this course