solvingnumericalanalysisproblemwithpython
Solving the Heat Equation Consider the heat equation ut = uxx, 0 ≤ x ≤ 1, t ≥ 0, with the initial condition u(0, x) = { 2x, 0 ≤ x ≤ 0.5 2 − 2x, 0.5 ≤ x ≤ 1 , (1) and boundary conditions u(t, 0) = 0, u(t, 1) = 0, t ≥ 0. A discretization of this system using forward differences in time and centered differences in space results in the recurrence uk+1i = uki + ∆t(∆x)2 ( uki+1 − 2uki + uki−1 ) , i = 1, · · · , n, where uki denotes the solution at spatial location i and time k. For more details, see Example 11.2 of the textbook. Using this full discretization with ∆x = 0.05 and ∆t = 0.0012, evolve the solution from t = 0 to t = 0.06. Plot the solution at the initial time, the final time, and periodically in between (say, every ten time steps or so). You should either use a 3d plot such as matplotlib’s plot_wireframe or plot the solutions for different times on the same set of axes and use different colored lines and a legend to distinguish the solution at different times. Repeat the previous part, but use the time step of ∆t = 0.0013. In a new figure, plot the solution in the same manner as before. Solve the same equation again, this time using the implicit method given in Section 11.2.2 of the textbook based on the backward Euler method. The recurrence for this method is uk+1i = uki + ∆t(∆x)2 ( uk+1i+1 − 2uk+1i + uk+1i−1 ) , i = 1, · · · , n. Again use ∆x = 0.05, but try a much larger time step of ∆t = 0.005 to advance the solution from t = 0 to t = 0.06. In a new figure, plot the solution at each time step. Solve the same equation again, this time using the Crank-Nicolson method given in Section 11.2.2. The recurrence for this method is uk+1i = uki + ∆t2(∆x)2 ( uk+1i+1 − 2uk+1i + uk+1i−1 + uki+1 − 2uki + uki−1 ) , i = 1, · · · , n. Again use ∆x = 0.05 and ∆t = 0.005 to advance the solution from t = 0 to t = 0.06. In a new figure, plot the solution at each time step. 1 Form the semidiscrete system for this equation given in Section 11.2.1 of the textbook. This discretization results in the following system of ODEs y′i(t) = 1(∆x)2 (yi+1(t) − 2yi(t) + yi−1(t)) , i = 1, · · · , n. Form the system with a finite difference spatial discretization with ∆x = 0.05 and use scipy.integrate.odeint to integrate it from t = 0 to t = 0.06 with ∆t = 0.005. In a new figure, plot the solution as before. Think about any differences in your plots and why they would occur given what you have learned about methods for solving PDEs. • Your code should produce five figures as described above. • Store the last step of every solution (i.e. the solution at t = 0.06) in the arrays solution_full_12, solution_full_13, solution_implicit, solution_cn, and solution_ode. 2 Solving the Heat Equation Solving the Wave Equation Use the method of lines and scipy.integrate.odeint to solve the wave equation utt = uxx, 0 ≤ x ≤ 1, t ≥ 0, with initial conditions u(0, x) = sin(πx), ut(0, x) = 0, 0 ≤ x ≤ 1, and boundary conditions u(t, 0) = 0, u(t, 1) = 0, t ≥ 0. Integrate from t = 0 to t = 2 for some reasonable value of ∆x. Plot the computed solution as a three-dimensional surface over the (t, x) plane. Note that the solution is periodic with period 2, so the solution should be the same at t = 0 and t = 2. You should note whether this is true for your computed solution, though you do not need to output anything from this check. Determine the maximum absolute error in the computed solution by comparing with the exact solution u(t, x) = cos(πt) sin(πx). Compute this error for solutions resulting from ∆x = 110·2∗∗k for k = 0, . . . , 5. On a log-log scale, plot the maximum error as a function of ∆x. Try to characterize the error as a function of ∆x. Store the errors in the array errors ordered by increasing values of k. • Plot the solution to the PDE with the given initial and boundary conditions for some reasonable value of ∆x. • Plot the maximum error as described above. Be sure to use a log-log scale. • As always, use appropriate titles, labels, etc. • Store errors as described above. 1 Solving the Wave Equation