代做MATH0033编程代写、代写Matlab程序语言、Matlab程序设计调试

日期：2021-01-04 11:27

MATH0033 Numerical Methods, 2020-2021, David Hewett

Computational homework 3

Differential equations

Exercises 1 and 2(a)-(e) (marked *) to be handed in and assessed

Deadline: 2200 GMT Sunday 10th January.

Please submit your solutions using the link on the course Moodle page.

You should submit a single pdf file, created using the Matlab publish

command, formatted as per the template provided on the Moodle page.

EXERCISE 1* Consider the initial value problem,

the exact solution of which is y(t) = tanh t = (e2t ? 1)/(e2t + 1).

(a) Compute numerical solutions to the initial value problem using the forward Euler method

with equally spaced points tn = n

20

N

, n = 0, . . . , N, for N = 19, 21, 40, 80 and 160.

Produce a plot comparing the numerical solutions with the exact solution.

For each choice of N compute the error maxn=0,...,N |un ? y(tn)| and store these errors in

an array efe. Also, compute and store the error |uN ? y(20)| at the end of the interval.

(b) Same as a) but using Heun’s method, to compute errors eheun.

(c) Produce a loglog plot of efe and eheun against the relevant N values and use your results

to read off the approximate convergence orders for the two methods. How do these

compare to the theory from lectures? (Hint: follow the approach we used in computational

homework 2)

(d) The exact solution y(t) has a horizontal asymptote y(t) → 1 as t → ∞. Does every approximation

obtained using the methods above reproduce the same asymptotic behaviour?

How well is the value of y(20) approximated?

Guided by the theory from lectures, but considering only the maximum value of |fy| on

the solution trajectory (rather than over all (t, y) ∈ R+ ×R), we might expect the critical

value of h below which the methods are stable to be h = 1. Is this what you observe in

your numerical results?

How does the lack of stability for N = 19 (h = 20/19 > 1) manifest itself for the two

different methods (forward Euler and Heun)?

1

EXERCISE 2 In this exercise we bring together some of the techniques we have studied in

the course to solve an initial boundary value problem (IBVP) involving a parabolic partial

differential equation, the heat equation.

The problem is: given T > 0 find u : [0, 1] × [0, T] 7→ R such that

?u

?t (x, t) ?

?

2u

?x2

(x, t) = 1 in [0, 1] × [0, T], (1)

with u(0, t) = u(1, t) = 0 for all t ∈ [0, T] and u(x, 0) = sin(πx) + 1

2

x(1 ? x). This problem is

a simple model for the evolution of the temperature distribution u in a metal bar undergoing

uniform heating, starting from an initial temperature distribution u(x, 0), where the bar is

held at constant (zero) temperature at its endpoints. For reference, the exact solution to this

problem is given by

u(x) = e

?π

2

t

sin(πx) + 1

2

x(1 ? x).

To solve the IBVP numerically we will use the so-called method of lines. This involves

discretizing the spatial differential operator ?

2u/?x2

in the same way for all times t, so that

the IBVP becomes a finite-dimensional system of ordinary differential equations in t, which

can then be solved using a time-stepping scheme such as forward Euler, backward Euler

or Crank-Nicolson. For the spatial discretization we shall use a standard finite difference

approximation. Given N > 0, define equally spaced nodes 0 = x0 < x1 < . . . < xN = 1 on

the interval [0, 1] by xn = nh, n = 0, . . . , N, where h = 1/N. For a function u(x, t) we denote

the approximation of u(xn, t) by un(t), and we replace the second derivative ?

2u/?x2

(x, t) by

the finite difference formula (cf. theoretical exercise sheet 1)

D2un(t) := un+1(t) ? 2un(t) + un?1(t)

h

2

.

This approximation turns the IBVP into an IVP involving a system of ordinary differential

equations satisfied by the functions un(t), n = 1, . . . , N ? 1. Explicitly, we obtain the system

du

dt (t) = f(t, u(t)) = ?Au + b, t ∈ [0, T], (2)

where u = (u1, . . . , uN?1)

T

, and, for a given N > 0, the matrix A and vector b can be

obtained in Matlab using the commands

>> h = 1./N;

>> A= (2/h^2)*diag(ones(N-1,1)) - (1/h^2)*diag(ones(N-2,1),1)...

- (1/h^2)*diag(ones(N-2,1),-1);

>> b = ones((N-1),1);

(a)* Consider the temporal discretization of the IVP system (2) using the forward Euler

method, the backward Euler method and Crank-Nicolson method, on a time mesh 0 =

t0 < t1 < . . . < tM = T, where M is the number of subintervals and tm = mk, m =

0, . . . , M, where k = T /M is the (uniform) temporal step size. (Note that I am using

k as the temporal step size and h as the spatial step size - don’t confuse the two!)

is the numerical approximation to u(xn, tm).

For each of the three methods write down the recurrence relation (in terms of the matrix

A) that has to be solved at each timestep to determine u

m+1 from u

m.

2

(b)* Now write a Matlab code which implements the three methods.

Hint: For the two implicit methods you need to solve a linear system at each timestep. For

simplicity I suggest you use Matlab’s backslash command, which calls a general-purpose

direct method. But you should be aware that in large-scale simulations one would use

either a special direct solver like the Thomas algorithm, or an iterative solver.

Note: I am suggesting that you implement the recurrence relations directly in your code,

rather than using the function files feuler.m etc. that I supplied. If you want to use the

latter, be warned that these codes only work for scalar problems, so you would need to

modify them appropriately to work for systems.

Run your code with the parameter choices N = 40 (giving a spatial step length h = 1/40),

k = h, and T = 0.2. Produce three plots (one for each method) showing the initial solution

and the solution at all of the time steps up to and including the final time T = 0.2 (on

the same axes, using the hold on command). Produce a further plot showing snapshots

of the exact solution at the same time steps.

What can be said qualitatively about the performance of the three methods, just by

studying your plots?

Hint: For forward Euler I suggest using the ylim command to adjust the vertical scale as

you see fit. The “vertical zoom” feature of the Matlab plot window is also helpful.

represent approximations to the exact solution u(x, t) on the INTERIOR spatial grid

points x1, . . . , xN?1.

N = 0 to each end of these vectors so you can plot on the whole spatial

interval [0, 1].

(c)* For each of the three methods, and for each N ∈ {10, 20, 40, 80, 160}, compute the rootmean-squared

error of the numerical approximation at the final time T = 0.2 using the

formula

E =

vuuth

N

X?1

n=1

(uMn ? u(xn, T))2,

where, as above, u

m

n denotes the approximation at time step m and in space node n, and

u(x, t) is the exact solution stated at the start of the question. As in part (b), use a time

step k equal to h.

Visualize your results on a loglog plot and, for the methods that converge, determine

the approximate order of convergence p for which E ≈ Chp

for some C > 0. (Hint: follow

the approach we used in computational homework 2)

(Bonus unassessed theoretical question: why do we include the scaling factor h inside the

square root in the definition of E?)

(d)*Now remove the assumption that h = k and, by varying h and k appropriately, investigate

the convergence order of the two implicit methods (backward Euler and Crank-Nicolson)

in the joint limit h → 0 and k → 0. That is, for each method determine constants p and

q for which the root-mean-squared error at T = 0.2 satisfies, for some C > 0, the bound

E ≤ C(h

p + k

q

).

Hint: to determine p, fix k and vary h, with k very small compared to all your values of h,

so that h

p dominates the bound. For instance, you could try M = 6400 and N = 10, 20, 40.

To determine q you would do the opposite: fix h very small compared to k (e.g. N = 640

and M = 8, 16, 32), so that k

q dominates the bound.

3

(e)* Now return to the forward Euler method. In lectures we proved that the matrix A is

symmetric positive definite, so that all its eigenvalues are real and positive. For N =

[10, 20, 40, 80, 160, 320] compute the eigenvalues of A using the eig command, and make

a precise conjecture about the behaviour of the maximum and minimum eigenvalues λmin

and λmax as the spatial step size h tends to zero.

Using your results, determine a condition on the temporal step size k of the form k ≤ Chp

for some C > 0 and p ≥ 1 (which you should specify), which ensures that the forward

Euler method is stable. (You will need to consider the spectral radius of the matrix

I ? kA.)

Verify that your stability criterion is correct by redoing your computations in (c) for the

forward Euler method, using your new stability criterion to choose k rather than setting

k = h. What convergence order do you observe as h → 0?

which can be approximated in Matlab using (provided N ? 1 is divisible by 3)

>> u0=[zeros((N-1)/3,1);ones((N-1)/3,1);zeros((N-1)/3,1)];

Let N = 22 and solve for one timestep, taking the timestep k equal to 1/N, using the

forward and backward Euler methods and the Crank-Nicolson method. Repeat using the

timesteps 1/(10N) and 1/(50N). Plot the approximations obtained and comment on the

results.

(g) (Bonus unassessed theoretical question!) At the start of the question I gave you an

analytical solution of the IVBP. Can you prove that this solution is the only solution? That

is, that the solution of (1) is unique? Hint: Try using the following “energy argument”.

Suppose that the IBVP has two solutions u1 and u2. Show that the difference v = u1 ? u2

satisfies the same IBVP, but with zero right-hand side and zero initial data. Using this

fact, along with integration by parts, show that the energy E(t) = (1/2) R 1

0

(v(t, x))2 dx

is a non-increasing function of t. Since E(0) = 0, deduce that E(t) = 0 for all t, and

conclude that v = 0, i.e. u1 = u2. You may assume that u1 and u2 are both continuous

functions of x.

4