Page - 177 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

5.2 Exercises 177 def K(u, t): N = len(u) - 1 K = zeros((N+1,N+1)) K[0,0] = 0 for i in range(1, N): K[i,i-1] = beta/dx**2 K[i,i] = -2*beta/dx**2 K[i,i+1] = beta/dx**2 K[N,N-1] = (beta/dx**2)*2 K[N,N] = (beta/dx**2)*(-2) return K import odespy solver = odespy.BackwardEuler(rhs, f_is_linear=True, jac=K) solver = odespy.ThetaRule(rhs, f_is_linear=True, jac=K, theta=0.5) solver.set_initial_condition(U_0) T = 1*60*60 N_t = int(round(T/float(dt))) time_points = linspace(0, T, N_t+1) u, t = solver.solve(time_points) Thefilerod_BE.pyhas all the details and shows amovieof the solution. Wecan run itwithany twewant, its size just impacts theaccuracyof thefirst steps. Odespysolversapplydensematrices! Looking at the entries of theKmatrix, we realize that there are at maximum three entries different from zero in each row. Therefore, most of the entries are zeroes. The Odespy solvers expect dense square matrices as input, here with .N C1/ .N C1/ elements. When solving the linear systems, a lot of storage andwork are spent on the zero entries in thematrix. Itwould bemuch moreefficient to store thematrixasa tridiagonalmatrixandapplya specialized Gaussian elimination solver for tridiagonal systems. Actually, this reduces the work fromtheorderN3 to theorderN . Inone-dimensionaldiffusionproblems, thesavingsofusingatridiagonalma- trixaremodest inpractice, since thematricesareverysmallanyway. Intwo-and three-dimensionalPDEproblems,however,onecannot afforddense squarema- trices. Rather, onemust resort tomore efficient storage formats and algorithms tailored to such formats,but this isbeyondthe scopeof thepresent text. 5.2 Exercises Exercise5.1: Simulateadiffusionequationbyhand Consider the problem given by (5.9), (5.10) and (5.14). SetN D 2 and com- puteu0i , u 1 i andu 2 i by hand for i D 0;1;2. Use these values to construct a test function for checking that the implementation is correct. Copy useful functions from test_diffusion_pde_exact_linear.py and make a new test function test_diffusion_hand_calculation. Filename:test_rod_hand_calculations.py.