COSC2500 Numerical Methods in Computational Science - Assignment 3: Linear DEs
MODULE3. ODES 63
Exercises—due 3:00 pm Thursday 19th October Required
The grade of 1–5 awarded for these exercises is based on the required exercises. If all of the required exercises are completed correctly, a grade of 5 will be obtained.
-
R3.1 From Sauer computer problems 6.1:
Using Euler’s method and a 4th-order Runge–Kutta method, solve the fol-lowing IVPs on the interval [0, 1] where y(0) = 1 and
a) y′ = t
b) y′=2(t+1)yc) y′ = 1/y2
for step sizes h = 0.1×2−k, for 0 ≤ k ≤ 5 (i.e., h = 0.1,0.05,0.025,...). Plot the solutions for h = 0.1, 0.05, 0.025, comparing with the exact solution. Plot a log–log plot of the error at t = 1 as a function of the step size h.
Also compare with the solution obtained using Matlab’s ode45 or another adaptive-step Runge–Kutta code.
A fixed-step Runge–Kutta code, RK4.m, is available on Blackboard.
-
R3.2 Consider the systems of linear DEs from Sauer computer problems 7.1–7.3, questions 1–2:
a)
b)
a) b)
i. y′′=y+32et,y(0)=0,y(1)=13e. ii. y′′ = (2+4t2)y, y(0) = 1, y(1) = e.
i. y′′ = 3y−2y′, y(0) = e3, y(1) = 1. Solve these DEs using the shooting method.
Solve these DEs using the finite difference method. Use the finite dif- ference substitutions for the derivatives to convert the DEs to linear systems, and solve the linear systems using the backslash operator or otherwise.
R3.3 Discuss and critique the following Matlab and C codes. Briefly compare them with the code for Matlab’s ode45 (you can display it in the command window using type ode45 or open it in the editor with open ode45).
function [t,y,last] = RK4(f,tspan,yi,dt)
%Standard RK4 algorithm.
%This code is a standard RK4 algorithm with a
%fixed step-size.
%It can be used to solve explicit first-order ODEs.
%The local truncation error is O(dt^5).
%"Time" is the name for the x-axis variable.
%
%Input:
% % % % %
f is a function handle for the first-order ODE
- dy/dt = f(t,y)
- f is assumed to be a row vector.
tspan is a vector that contains ti and tf, e.g. tspan = [ti,tf]
64
COSC2500/COSC7500—Semester 2, 2023
ti is the initial time
tf is the final time
yi is the initial values for the ODE.
- yi is assumed to be a vector. dt is the step size
%
%
%
%
%
%
%Output:
%
%
%
%
%
%
%Hint - call the function like this:
% [T,Y] = RK4(f,tspan,yi,dt)
t is a vector of time values.
y is the approximate solution curve for the
initial value problem of the ODE. last is the last y values - useful if you are
using the last values for something
ti = tspan(1); tf = tspan(2);
num_steps = ceil((tf-ti)/dt);
%ceil forces it to be an integer
t = linspace(ti,tf,num_steps+1).’;
%creates time vector, then transposes
if size(yi,2) > 1 yi = yi.’;
end
%Initialising y
y = [yi,zeros(size(yi,1),num_steps)];
%Application of the RK4 algorithm.
for
end
y =
last = y(end,:);
end
void runge4(double x, double y[], double step) {
n = 1:num_steps
k1 = dt*f(t(n),y(:,n));
k2 = dt*f(t(n) + 0.5*dt, y(:,n) + 0.5*k1);
k3 = dt*f(t(n) + 0.5*dt, y(:,n) + 0.5*k2);
k4 = dt*f(t(n) + dt, y(:,n) + k3);
y(:,n+1) = y(:,n) + (1/6)*(k1+ 2*k2 + 2*k3 + k4);
y.’;
double h=step/2.0,
t1[N], t2[N], t3[N],
k1[N], k2[N], k3[N],k4[N];
int i;
/* the midpoint */
for (i=0;i<N;i++) {t1[i]=y[i]+0.5*(k1[i]=step*f(x, y, i));}
MODULE3. ODES 65
for (i=0;i<N;i++) {t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));}
for (i=0;i<N;i++)
{t3[i]=y[i]+ (k3[i]=step*f(x+h, t2, i));}
for (i=0;i<N;i++)
{k4[i]= step*f(x+step, t3, i);}
for (i=0;i<N;i++) {y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;}
}
double f(double x, double y[], int i) {
/* derivative of first equation */
if (i==0) return(exp(-2*x) - 3*y[0]);
if (i!=0) exit(2);
}
-
R3.4 Discuss, with critical analysis, a paper from the research literature.
A list of papers will be available on Blackboard if you have difficulty findinga suitable paper to discuss.
-
R3.5 Find a median number from a unsorted array in C language.
For example, a sequence is 12, 10, 15, 3, 23, the median will be 12 because it is the middle element of the sorted sequence 3, 10, 12, 15, 23. A sequence is 12,10, 15, 3, the median will be 11 because it is the middle element of the sorted sequence 3, 10, 12, 15 since (10 + 12)/2 = 11.
-
a) The size and the elements of the array should be input from console. The first message asks about the size of the array and the second one asks about the elements of the array. For example, if the array has 4 elements as: 12,10, 15, 3. The console should be:
Enter the size of the array: 4 Enter the elements of the array: 12 10 15 3 -
b) Write a function to find the median numbers using sort algorithms. The sort algorithm can be found from any textbooks and web re- sources, such as
i. Selection Sort (With Code in Python/C++/Java/C) https://www. programiz.com/dsa/selection-sort
ii. C Program: Insertion sort algorithm https://www.w3resource. com/c-programming-exercises/searching-and-sorting/ c-search-and-sorting-exercise-4.php
-
c) Write a function to find the median number without sorting the se- quence (such as finding the k largest numbers or other methods).
-
d) Test the above two functions in a main function to verify that both functions are correct. Show that your code works for both even and odd size arrays.
-
66 COSC2500/COSC7500—Semester 2, 2023
Additional
Attempts at these exercises can earn additional marks, but will not count towards the grade of 1–5 for the exercises. Completing all of these exercises does not mean that 6 marks will be obtained—the marks depend on the quality of the answers. It is possible to earn all 6 marks without completing all of these additional exercises.
-
A3.1 Estimate the error in a numerical solution obtained using Euler’s method from the convergence of the solution as the step size is changed. Then com- pare your estimate of the error with actual error determined from compar- ison with the analytical solution.
-
A3.2 Below is a model that is a variation of the Lotka–Volterra equations (i.e., predator–prey equations):
-
x′ = 1.3x − 0.5xy
-
y′ = 0.2xy − 0.5y
Choose an ODE solver and solve this system for both x0 and y0 equal to 12 betweent=0andt=500. Plotxandyvstandthenax-yplot. Briefly comment on the behaviour of the system and whether this model seems realistic (e.g., to model a system consisting of predators and prey).
One possible way to solve a system in Matlab using ode45 is to create an anonymous function that returns a column vector. To make it work with most of the Matlab ODE solvers, the inputs should be t and an arbitrary value (say) u. Every incidence of x is replaced with the indexed variable u(1) and every incidence of y is replaced with u(2).
For reference on how you would represent the system in Matlab, the fol- lowing example is provided. For the system
x′ = xy+x y′ = y
the corresponding anonymous function code is
func = @(t,u) [u(1)*u(2) + u(1); u(2)];
% You should not be solving this equation % for the module exercise,
% this is an example... -
-
A3.3 Investigate Sauer computer simulation 6.3.2 (the pendulum). For some spe- cific tasks, see computer problems 6.3.
-
A3.4 Investigate Sauer computer simulation 6.3.3 (orbital mechanics). For some specific tasks, see computer problems 6.3.
-
A3.5 Sauerrealitycheck6—TheTacomaNarrowsbridge(Sauerpg322in2E,328 in 1E). As well as the suggested activities, try using Matlab’s ode45 solver.
-
A3.6 Simulate the motion of a paper airplane.
(Googling for "paper airplane" "differential equations" lift drag will provide a convenient starting point.)
MODULE3. ODES 67
A3.7 Considerthesystemsofnon-linearDEsfromSauercomputerproblems7.1– 7.3, questions 3–4:
a) i. y′′=18y2,y(1)=1,y(2)= 1. ′′ −2y 23 12
ii. y =2e (1−t ),y(0)=0,y(1)=ln2. b) i. y′′ =ey,y(0)=1,y(1)=3.
ii. y′′ =siny′,y(0)=1,y(1)=−1.
Solve these using the finite difference method, using the non-linear code from Sauer or otherwise. You might like to compare your solutions with those obtained using the shooting method.
-
A3.3 Solve the DEs from the required exercises and/or A3.7 using the finite ele- ment method, collocation method, or similar method.
-
A3.4 Estimate the error in a numerical solution obtained using the finite dif- ference method from the convergence of the solution as the step size is changed. Then compare your estimate of the error with actual error de- termined from comparison with the analytical solution.
-
A3.5 Solve the system of equations in Gramotnev and Nieminen 1999. The ap- pendix in the paper describes a finite difference solution of the system. Can you solve it using the shooting method?
-
A3.6 Model a multi-species predator–prey relationship using a homogeneous linear matrix DE. If we begin with known populations of each species, we have an initial value problem; if If we know a mixture of final and initial populations, we have a BVP. (You might like to try N = 5.)
Programming hints
R3.2(a) For the shooting method, you’re combining and initial value problem solver and a root-finder. You can use whichever IVP solver (from Module 4 or elsewhere) and root-finder (from Module 2 or elsewhere) you prefer. The only tricks are that
-
Your IVP solver will want a system of 1st order ODEs. You are starting with a second order ODE. So, the first step is to convert the 2nd order ODE into a system of 2 1st order ODEs, as described at the start of Module 4.
-
You need to write a function that the root-finder will find the roots of. For the shooting method, we guess the missing initial condition, and solve the system of ODEs as an IVP. We can then compare our solution to the other boundary condition (which we haven’t used yet). We want the difference between our solutions and this boundary condition to be zero. That’s our root-finding problem. So we need a function like
function final_error = bc_mismatch(x)
% Find the roots of this for R5.1 shooting method
% Solve the IVP
[t,y]=ode45(@r51_1a,[0 1],[0 x]);
68
COSC2500/COSC7500—Semester 2, 2023
% Compare with final BC
bc = (1/3) * exp(1); final_error = y(end,1) - bc;
return
The parameter passed to this function, x, is your guess for the missing initial condition. Since your root-finder will only want to pass one parameter to the function it’s finding roots of, it’s easiest to hard-code things like the ODE system function name, the boundary conditions, etc. in the file as in this example. Depending on which order you write your 2 1st order ODEs in, you might need to swap the order of the 2 boundary conditions passed to your IVP solver, and to compare the final value y(end,2).
As far as your root-finder cares, this function is just another function to find roots of. You can also calculate values for different x and plot a graph; this will be easiest to do using a loop to cycle through the different values of x.
R3.2(b) If you want to solve these as linear systems, you’ll need to write your own code. Once you have your linear system, solving it is easy. There are two steps you will need to do first:
-
Substitute the finite difference approximations into your ODE. Here is an example, with a complicated ODE:
OurODE:y′′ =y′+tsint
Our finite difference approximation for y′′:y′′=(y −2y +y )/h2 n n+1 n n−1
Note that yn = y(tn) and yn+1 = y(tn+1) = y(tn + h).
Let’s use the forward difference for y′: y′n = (yn+1 − yn)/h. Finally, we’ll replace t by tn.
So, our ODE becomes:yn+1 −2yn +yn−1 = yn+1 −yn +tn sintn. h2 h
We want to write this in our usual form for a linear equation in y1, y2, etc., as An yn−1 + Bn yn + Cn yn+1 = Dn. Re-arranging, we get
1 yn−1 + 1 − 2 yn + 2 − 1 yn+1 = tn sin tn. h2 h h2 h2 h
-
Generate the linear system. The easy way is to put our initial and final boundary conditions as the 1st and last rows.
% How many points do we want?
N = 10;
% Pre-allocate memory for matrix and vector
A = zeros(N,N); c = zeros(N,1);
MODULE3. ODES 69
tn = linspace(0,1,N); % Need to change for 2(a), % since it goes to 1.5
h = t(2) - t(1);
% Boundary conditions
A(1,1) = 1;
c(1) = 0; % First BC
A(N,N) = 1;
c(N) = exp(1)/3; % End BC
% Make the matrix
for n = 2:(N-1)
A(n,n-1) = 1/h^2;
A(n,n) = 1/h - 2/h^2
A(n,n+1) = 2/h^2 - 1/h;
c(n) = tn(n) * sin(tn(n));
end
Finally, solve and plot!
