Questions:
---------1. (2 points) Develop a Runge-Kutta ODE solver function
named odex45x() as an updated version of the ode02()
solver of the lecture notes, or of the corresponding
solver in lecture transcripts. Your solver will use
the Runge-Kutta formula (with dx1, dx2, dx3 and dx4)
given in lecture L03_04_ODE_solvers_revisited.pdf.
Develop a script that uses your odex45x() to solve
the predator-prey problem discussed in lecture:
dH/dt = 0.5 H - 0.02 H P
dP/dt = - 0.4 P + 0.001 H P
for 100 years, with H(0) = 600 and P(0) = 10. Plot
the populations of H and P versus time, and the
H-P phase portrait on a single figure with 4
subplots where the upper row displays results
of your oderx45x() and the lower row displays results
obtained with MATLAB's ode45().
Paste your odex45x, the script that you used to solve
this problem and the figure that you produced in the
space below. Also, state the degree to which the
solution produced by your odex45x is the same (or not)
as that produced by ode45.
2. (2) Symplectic integrators are commonly used to
solve Hamiltonian systems for the position of bodies
evolving over time, under a force field. They have
applications in molecular dynamics (eg, Lennard-Jones
potential) and celestial mechanics (eg. gravitational
potential), among others. Much like Runge-Kutta and
other solvers, they advance a solution by one time
step based on solution values at a previous time
step. However, in distinction to conventional solvers,
they use different formulas for each state variable.
As they are specialized to calculations of motion,
they consider position (X), velocity (V) and acceleration
(A) as the state variables to be computed.
The Verlet method is an example of such an integrator.
It advances the solution from one time step to the
next using the formulas:
X(t+dt) = X(t) + V(t)*dt + 0.5*A(t)*dt*dt eq. 2a
A(t+dt) = F(t+dt) / M
eq. 2b
V(t+dt) = V(t) + 0.5*(A(t+dt)+A(t))*dt
eq. 2c
where F(t) is the force (or sum of all forces),
acting at time t, on the body whose motion is
being computed, and M is the mass of this body.
The state variables, X, A, V, are vectors (with 2
components in 2-D or 3 components in 3-D), and the
force may be obtained from an appropriate function
that returns a vector (with 2 items in 2-D, eg. Fx
and Fy, or 3 items in 3-D).
From a computational perspective, if a single body
is in motion, then each of its state variables can
be stored in a matrix with as many rows as there
are times at which position is to be computed, and
as many columns as there are spatial dimensions in
the problem being solved, for example (in 2-D with
n time steps):
[ x1 y1] [vx1 vy1] [ax1 ay1]
[ x2 y2] [vx2 vy2] [ax2 ay2]
X = [ x3 y3], V = [vx3 vy3], A = [ax3 ay3]
[... ...] [... ...] [... ...]
[ xn yn] [vxn vyn] [axn ayn]
Develop a simplectic solver named hsi02() that
takes a force function, time vector, body mass,
initial position vector and initial velocity vector
as inputs and returns the time vector (as a column
vector), and the matrices of positions, velocities
and accelerations computed for a single body moving
under the effect of the specified force. Your
hsi02() should work for both 2-D and 3-D problems.
A possible skeleton is as follows:
function [t,X,V,A] = hsi02(force,tv,M,X0,V0)
t = tv(:);% output time vector (column)
nd = ___;
% number of spatial dimensions
nt = ___;
% number of time steps
X = zeros(__,__); % position matrix (for solution)
V = zeros(__,__); % velocity matrix (for solution)
A = zeros(__,__); % acceleration matrix (for solution)
X(1,:) = ___;
% store initial position in sol matrix
V(1,:) = ___;
% store initial velocity in sol matrix
A(1,:) = force(X(__))/M;
% store initial accel in sol
for k = ______
X(...) = ...;
% update position for time step k
A(...) = ...;
% update acceleration for time step k
V(...) = ...;
% update velocity for time step k
end
Test your hsi02() using the following script:
% -------- start of test script --------------------% problem constants
G = 6.674e-11; % gravitational constant (N m^2/Kg^2)
Me = 5.97e24; % mass of earth (Kg)
Mm = 7.35e22; % mass of moon (Kg)
Dem = 378000e3; % earth-moon distance (m)
Vm = 1.02e3; % moon's orbital speed (m/s)
Fgra = @(X) -(G*Me*Mm./sum(X.^2))*X/sqrt(sum(X.^2)); % gravity
% 2-D case: set initial values and time vector
X0 = [Dem 0]; % initial moon position ( x, y) (m)
V0 = [ 0 Vm]; % initial moon velocity (vx,vy) (m/s)
tv = 24*3600*linspace(0,120,241); % 120 days (in seconds)
% 2-D case: solve and plot
[t, Xm] = hsi02(Fgra,tv,Mm,X0,V0);
figure(1); plot(Xm(:,1)/1000,Xm(:,2)/1000)
xlabel('x (km)'); ylabel('y (km)');
axis equal; axis([-1 1 -1 1]*1.1*Dem/1000)
% 3-D case: set initial values and time vector
X0 = [Dem 0 0]; % initial moon position (m)
V0 = [ 0 Vm 0]; % initial moon velocity (m/s)
% 3-D case: solve and plot
[t, Xm] = hsi02(Fgra,tv,Mm,X0,V0);
figure(2); plot3(Xm(:,1)/1000,Xm(:,2)/1000,Xm(:,3)/1000)
xlabel('x (km)'); ylabel('y (km)'); zlabel('z (km)');
axis equal; axis([-1 1 -1 1 -1 1]*1.1*Dem/1000)
% -------- end of script ----------------------Paste your hsi02() and the two figures produced by
the test script in the space below.
3. (2) An approximate French flag (vertical tricolour,
blue, white, red (equally spaced), with height:width ratio
of 2:3, can be constructed in MATLAB using a 3-D array of
uint8 data items (vertical-position as rows, horizontal-position
as columns, color in 3 planes: Red, Green and Blue intensities,
each 0 to 255), and then displayed, using the script:
Fflag = uint8(255*ones(100,150,3));
Fflag(:, 1:50, 1:2) = 0;
Fflag(:,101:150,2:3) = 0;
imshow(Fflag)
Develop a similar script to produce an approximate flag
for Benin, with 2:3 height:width ratio, a vertical green
band on the left third and horizontal yellow and red
bands on the right third. Paste your script and flag in
the space below.
4. (2) Consider the concentration field resulting from
solute transport under advective and diffusive processes
in three spatial dimensions (see Question 5 of HW #05
and HW #06 and Question 4 of HW #08 for the 2-D case):
c(x,y,z,t) = [M/(4 pi D t)^(3/2)]
* e^{-[(x - vx t)^2 + y^2 + (z - vz t)^2]/(4 D t)}
with M = 10 grams, vx = 0.5 cm/min, vz = 0.3 cm/min
and D = 0.25 cm^2/min.
Use MATLAB to plot slices of this concentration field
at time t = 0.5, t = 1, t = 2 and t = 5 minutes,
respectively (as in prior HW), using a for-loop
and 4 subplots (1 subplot per time). Use x-y-z grids
where x goes from -0.5 cm to 3.5 cm, y ranges from
-1.0 to 1.0 cm and z goes from -0.5 to 2.5 cm. Place
a slice at z = 0, a slice at y = 0, a slice at x = 0
and another slice at x = 2.5 cm. (be sure to
use flat shading and equal x-y-z axis scales). Label
your axes, put a title and a colorbar on each plot.
Paste the script that you developed for this task
in the space below and also paste the figure you
produced in this space.
5. (2) Use MATLAB to produce level surfaces (plots) of the
3-D diffusive field of the previous question, with the
same parameter values, times and axes, except for the
y-axis that should now go from 0.0 to 1.0 cm. Use five
or more concentration levels such that the nature of the
concentration field is readily ascertained from your
plots. Paste your script and resulting figure in the
space below.
Purchase answer to see full
attachment