写作MATH0033编程、辅导Matlab

” 写作MATH0033编程、辅导MatlabMATH0033 Numerical Methods, 2020-2021, David HewettComputational homework 3Differential equationsExercises 1 and 2(a)-(e) (marked *) to be handed in and assessedDeadline: 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 publishcommand, 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 methodwith equally spaced points tn = n20N, 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 inan array efe. Also, compute and store the error |uN y(20)| at the end of the interval.(b) Same as a) but using Heuns method, to compute errors eheun.(c) Produce a loglog plot of efe and eheun against the relevant N values and use your resultsto read off the approximate convergence orders for the two methods. How do thesecompare to the theory from lectures? (Hint: follow the approach we used in computationalhomework 2)(d) The exact solution y(t) has a horizontal asymptote y(t) 1 as t . Does every approximationobtained 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| onthe solution trajectory (rather than over all (t, y) R+ R), we might expect the criticalvalue of h Below which the methods are stable to be h = 1. Is this what you observe inyour numerical results?How does the lack of stability for N = 19 (h = 20/19 1) manifest itself for the twodifferent methods (forward Euler and Heun)?1EXERCISE 2 In this exercise we bring together some of the techniques we have studied inthe course to solve an initial boundary value problem (IBVP) involving a parabolic partialdifferential equation, the heat equation.The problem is: given T 0 find u : [0, 1] [0, T] 7 R such thatut (x, t) 2ux2(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) + 12x(1 x). This problem isa simple model for the evolution of the temperature distribution u in a metal bar undergoinguniform heating, starting from an initial temperature distribution u(x, 0), where the bar isheld at constant (zero) temperature at its endpoints. For reference, the exact solution to thisproblem is given byu(x) = e2tsin(x) + 12x(1 x).To Solve the IBVP numerically we will use the so-called method of lines. This involvesdiscretizing the spatial differential operator 2u/x2in the same way for all times t, so thatthe IBVP becomes a finite-dimensional system of ordinary differential equations in t, whichcan then be solved using a time-stepping scheme such as forward Euler, backward Euleror Crank-Nicolson. For the spatial discretization we shall use a standard finite differenceapproximation. Given N 0, define equally spaced nodes 0 = x0 x1 . . . xN = 1 onthe interval [0, 1] by xn = nh, n = 0, . . . , N, where h = 1/N. For a function u(x, t) we denotethe approximation of u(xn, t) by un(t), and we replace the second derivative 2u/x2(x, t) bythe finite difference formula (cf. theoretical exercise sheet 1)D2un(t) := un+1(t) 2un(t) + un1(t)h2.This approximation Turns the IBVP into an IVP involving a system of ordinary differentialequations satisfied by the functions un(t), n = 1, . . . , N 1. Explicitly, we obtain the systemdudt (t) = f(t, u(t)) = Au + b, t [0, T], (2)where u = (u1, . . . , uN1)T, and, for a given N 0, the matrix A and vector b can beobtained 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 Eulermethod, 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 usingk as the temporal step size and h as the spatial step size – dont 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 matrixA) that has to be solved at each timestep to determine um+1 from um.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. Forsimplicity I suggest you use Matlabs backslash command, which calls a general-purposedirect method. But you should be aware that in large-scale simulations one would useeither 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 thelatter, be warned that these codes only work for scalar problems, so you would need tomodify 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 solutionand the solution at all Of the time steps up to and including the final time T = 0.2 (onthe same axes, using the hold on command). Produce a further plot showing snapshotsof the exact solution at the same time steps.What can be said Qualitatively about the performance of the three methods, just bystudying your plots?Hint: For forward Euler I suggest using the ylim command to adjust the vertical scale asyou 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 gridpoints x1, . . . , xN1.N = 0 to each End of these vectors so you can plot on the whole spatialinterval [0, 1].(c)* For each of the three methods, and for each N {10, 20, 40, 80, 160}, compute the rootmean-squarederror of the numerical approximation at the final time T = 0.2 using theformulaE =vuuthNX1n=1(uMn u(xn, T))2,where, as above, umn denotes the approximation at time step m and in space node n, andu(x, t) is the exact solution stated at the start of the question. As in part (b), use a timestep k equal to h.Visualize your results on a loglog plot and, for the methods that converge, determinethe approximate order of convergence p for which E Chpfor some C 0. (Hint: followthe approach we used in computational homework 2)(Bonus unassessed theoretical question: why do we include the scaling factor h inside thesquare root in the definition of E?)(d)*Now remove the assumption that h = k and, by varying h and k appropriately, investigatethe 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 andq for which the root-mean-squared error at T = 0.2 satisfies, for some C 0, the boundE C(hp + kq).Hint: to Determine p, fix k and vary h, with k very small compared to all your values of h,so that hp 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 = 640and M = 8, 16, 32), so that kq dominates the bound.3(e)* Now return to the forward Euler method. In lectures we proved that the matrix A issymmetric 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 makea precise conjecture about the behaviour of the maximum and minimum eigenvalues minand 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 Chpfor some C 0 and p 1 (which you should specify), which ensures that the forwardEuler method is stable. (You will need to consider the spectral radius of the matrixI kA.)Verify that your stability criterion is correct by redoing your computations in (c) for theforward Euler method, using your new stability criterion to choose k rather than settingk = 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 theforward and backward Euler methods and the Crank-Nicolson method. Repeat using thetimesteps 1/(10N) and 1/(50N). Plot the approximations obtained and comment on theresults.(g) (Bonus Unassessed theoretical question!) At the start of the question I gave you ananalytical solution of the IVBP. Can you prove that this solution is the only solution? Thatis, 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 u2satisfies the same IBVP, but with zero right-hand side and zero initial data. Using thisfact, along With integration by parts, show that the energy E(t) = (1/2) R 10(v(t, x))2 dxis a non-increasing function of t. Since E(0) = 0, deduce that E(t) = 0 for all t, andconclude that v = 0, i.e. u1 = u2. You may assume that u1 and u2 are both continuousfunctions of x.如有需要，请加QQ：99515681 或WX：codehelp

“