本次cs代写的主要内容是MATLAB数学关于数值分析的
Due Date: Friday 26 March, 2021.
Total points: 49.
- Assignments below are either theory or computational. For the computational assignments you may use any programming language you want except Maple or Mathematica. Besides your written/typed solutions you must also submit your code.
- This assignment should be submitted to Crowdmark. Please make sure your upload has the correct orientation.
- Office hrs: Tuesday’s and Thursday’s 10:00–11:00am (local time Waterloo, ON). A link will be emailed to you before each office hour.
- If you want to use your one time 3-day extension, the due date is Wednesday 31 March, 2021. Please email the TA if you plan to use your extension.
Excersises X.Y are from the book “A First Couse in the Numerical Analysis of Differential Equations” by Arieh Iserles, 2nd edition.
- (0 points) Please sign the Academic Integrity Checklist. If you do not sign the Academic Integrity Checklist you will receive a 0 for this assignment.
- (a) (5 points) Exercise 3.8.
(b) (3 points) Explain whether the following three-stage IRK is or is not a collocation method:k1=yn+h1f(tn+1h,k1)+ 1f(tn+1h,k2)+ 1f(tn+3h,k3) 8 4 16 2 16 4
k2=yn+h3f(tn+1h,k1)+ 3f(tn+1h,k2)+ 7f(tn+3h,k3) 16 4 32 2 32 4
k3 =yn+h12f(tn+41h,k1)+18f(tn+12h,k2)+18f(tn+34h,k3) yn+1=yn+h3f(tn+1h,k1)+3f(tn+1h,k2)+ 1f(tn+3h,k3)
- (5 points) Exercise 4.4.
- (5 points) Exercise 4.6a. You do not need to do parts (b) and (c).
- To study the trajectories of the stars in the galaxy, H ́enon and Heiles developed the following system of equations
d2x = −∂V , d2y = −∂V , (1)
dt2 ∂x dt2 ∂y
where V (x, y) is the H ́enon-Heiles potential given by
V (x, y) = 1 x2 + y2 + x2y − 31 y3 .
2 dt dt
this energy is conserved during motion.
(a) (1 point) Show that the H ́enon-Heiles system (1) can be written as the following first order system:
10 4 5 2 10 4
2 TheenergyofthesystemisdefinedasE(t)=V(x(t),y(t))+1 (dx)2 +(dy)2. Itcanbeshownthat
dx =p dt
dp =−x−2xy dt
dy =q dt
dq =−y−(x2−y2)
(b) (7 points) Use any 4th order explicit numerical method to find an approximate solution to the H ́enon-Heiles problem (write down which method you use). Consider the following two sets of initial conditions:
and
x(0) = 0, x(0) = 0,
dx (0) = 0.436630946376151, dt
dx (0) = 0.427001854016272, dt
y(0) = 0.095,
y(0) = 0.095,
dy (0) = 0.03, dt
dy (0) = 0.096. dt
For each set of initial conditions verify that energy is indeed conserved (up to approximately nine digits behind the comma) on the time interval t ∈ [0,1000] (as time step, use h = 10−2). Also, plot the trajectory in the (x, y) plane for t ∈ [0, 1000].
6. In this question we will solve a model of flame propagation. This problem is described on Cleve Moler’s blog (Cleve Moler is the inventor of Matlab) and is restated here. If you light a match, the ball of flame grows rapidly until it reaches a critical size. Then it remains at that size because the amount of oxygen being consumed by the combustion in the interior of the ball balances the amount available through the surface. The model is given by the following non-linear ODE:
y′=y2−y3 y(0)=δ t∈[0,TF], TF =2δ,
where δ describes the initial radius of a ball and y(t) represents the radius of the ball. For very small
δ, this is a stiff ODE. For larger δ, the problem is not very stiff.
- (a) (7 points) Implement (explicit) Euler’s method for this problem. Take δ = 0.02. Plot the
solution when h = TF /50 and when h = TF /100. Explain the difference in solution.
- (b) (3 points) Take now δ = 0.0001. Plot the solution when h = TF /10000 and h = TF /20000.
Explain what you see.
From (b) we just saw that we need on the order of 20000 time steps for Euler’s method to be stable! It would be better to use an A-stable method. Below we will implement backward Euler’s method. Backward Euler’s method for the flame propagation model is given by
y = y + hy2 − hy3 , n = 0, 1, … (3) n+1 n n+1 n+1
We immediately see that we run into a problem here: this is a non-linear equation for yn+1. To solve this non-linear equation, we use Newton-Raphson.
(c) (3 points) Show that Newton-Raphson applied to (3) results in the following iterative scheme: w[i+1] = w[i] − F(w[i]) , i = 0,1,…, (4)
F′(w[i])
whereF(w)=w−yn−hw2+hw3 andF′(w)=1−2hw+3hw2. We implement Newton-Raphson as follows:
- First we need a good initial guess. Let’s choose w[0] = yn.
- We then start iterating over i in (4). We continue to iterate until |w[i+1] − w[i]| < 10−10.
- Once |w[i+1] − w[i]| < 10−10 is satisfied, we set yn+1 = w[i+1], and therefore we have solved (3) up to a tolerance of 10−10, which is usually “good enough”.
Pseudo code for Newton’s method applied to the Backward Euler method for the flame propagation problem, (3), is given by:
δ = (to be set by user) TF =2/δ
h = (to be set by user) n=1
t=0
yn = δ
while t < TF
Set w[0] = yn Set ε = 1
while ε > 10−10
Compute F (w[i]) = w[i] − yn − h(w[i])2 + h(w[i])3 Compute F ′(w[i]) = 1 − 2hw[i] + 3h(w[i])2 Compute w[i+1] = w[i] − F (w[i])/F ′(w[i]) Compute difference ε = |w[i+1] − w[i]|
Update w[i] = w[i+1] end
Update yn+1 = w[i+1] Update t = t + h Update n = n + 1
end
(d) (10 points) Implement backward Euler’s method as just described. Take δ = 0.0001 and plot the solutionforh=TF/20andh=TF/100. Giventhatweexpectthetransitionfromy=δtoy=1 to occur around t = 1/δ, we see that the solutions obtained with h = TF /20 and h = TF /100 are not very accurate. How small do we need to choose h such that we obtain, at least visually, a sensible solution?
We learn from this exercise that for stiff ODE’s it may be worth spending a bit more time implementing an implicit numerical scheme than trying to solve the problem with explicit numerical schemes. We saw from question (b) that we need on the order of 20000 time steps for Euler’s method to be stable. For backward Euler’s method we saw from question (c), that even at a time step of h = TF /20 the method is still stable. Not stability, but accuracy is now the determining factor in choosing our time step.
