Trinh @ Bath

A primer in numerical methods

Beyond visualizing the mathematical concepts in pattern formation, you may be curious to learn about how the animations were actually implemented using MATLAB. In this section I want introduce some of the key ideas behind the code that was used to generate the animated gifs, guide you through some of the code, and make it accessible for you to make similar animations as well.

First, I will go through the (explicit and implicit) numerical solution of the 1D Heat Equation, and then I will step through the explicit numerical solution for a 2D Heat Equation.

Note: I assume that you are familiar with the Forward and Backward Euler Methods to solve ordinary differential equations.

The Heat Equation

Recall that the Heat equation (also the diffusion equation) is: utDu2ux2=0,utDu2ux2=0, where uu is temperature, and DuDu is a constant that determines how fast heat spreads out. Before we go into the numerics let us step through the solution to the heat equation.

Imagine we have an insulated rod with initial condition u(x,0)=u0 boundary conditions u(x,t)=u(L,t)=0. We want to solve the PDE in Equation [???]. The trick is to assume that the solutions are separable, so u(x,t)=T(t)X(x)

We can substitute this solution into Equation [???] to get T(t)X(x)=ut=DuT(t)X(x)=2ux2.

Rearranging terms we get

T(t)DuT(t)=X(x)X(x).

Note that the LHS (left hand side) is a function of time where the RHS (right hand side) is a function of position. In order for this equality to be satisfied for all x and t, we must have that they are both equal to a constant. We let

T(t)DuT(t)=X(x)X(x)=λ2λ>0.

Solving for each variable we get that T(t)=Aeλ2DutX(x)=Bcos(λx)+Csin(λx)

Here, A is a positive real constant. We can apply our boundary conditions to find what B and C should be. Since the edges at x=0 is clamped at zero (u(0,t)=0), then B=0. The rod is also clamped at zero at x=L, and since C0, then sin(λL)=0 This means that λn=nπL where nZ.

Note that:

  • n=0 give a trivial solution
  • nZ just makes all the coefficients negative, so we can just consider the positive values of n.

As a result, we get the following eigenfunctions:

un(x,t)=Aneλ2nDutsin(λnx)n=1,2,..., with eigenvalues

λn=nπL

We are not done yet, because we still need to find the conditions s.t. we satisfy the initial conditions. Since linear combinations of solutions will still be solutions, then u(x,t)=n=1Aneλ2nDutsin(λnx), is still a solution. We would like it then, if we could express u0 as a sum of sine functions as shown in Equation [???]. This is true for most u0 and it is the theory of Fourier Series.

Although I will not get into the details here, it can be shown that the Fourier sine coefficients are:

An=2LL0u0(x)sin(λnx)dx

However, in general, it is not always easy to obtain solutions to PDE's in closed form, and arguably even Equation (???) is not as easy to understand (unless you can visualize an infinite sum!). Instead, we turn to numerical simulations to Partial Differential equations. One of the methods used is known as the finite-difference method. (For a more complete discussion and more Matlab examples and exercises see [(kharab_2002)]).

Forward (explicit) method

One of the difficulties when numerically simulating the heat equation is that (???) contains partial derivative in both time and space. We can approximate the partial derivatives of u(x,t) in Eqn. (???) as follows: 2ux2u(x+Δx,t)2u(x,t)+u(xΔx,t)Δx2,utu(x,t+Δt)u(x,t)Δt.

This finite difference method is also known as the Forward-Time Central-Space method, because at each iteration we approximate the value of u by considering considering the change in space centered around a point in space x using values at x±Δx, and the change in time as in the forward Euler method. If you were to do it manually you would have a “T-shaped stencil” that you could move accross the discretized 1D rod to find the values at the next time step as shown in the figure below.

 Forward finite difference stencil

Let n be the index for the time steps, and m be the index for position. We'll also use k=Δt and h=Δx, and write u(m,n)=unm

Then, we can solve the heat equation numerically by using: 1k(un+1munm)=1h2(unm12unm+unm+1)un+1m=unm+kh2(unm12unm+unm+1)un+1m=kh2unm1+(12kh2)unm+kh2unm+1.

What happens at the edges?

Let us discretize our rod at M points, un1,un2,...,unM1,unM. Notice that stencil will work at every position along a one-dimensional rod except for the two edges (at u1 and uM), because we have no a priori values for u0 and uM+1 to approximate the centered finite difference. We can work around this by using the fact that we applied periodic boundary conditions. This means that un1=unM, so: un0=unM1,unM+1=un2

If you use MATLAB to program the numerical solution to the 1D Heat Equation, it is convenient to rewrite (???) as a matrix operation, un+1=Bun, where un are the values of u at time n for all values x1,...,xM, and B=(ab0b0bab000ba00bab0b0ba), where a=12k2, and b=kh2, and where B is a tridiagonal matrix except at B1,M1=BM,2=b that satisfies the boundary conditions.

Backward (implicit) method

The forward method works well, but is not guaranteed stability. Implicit methods are guaranteed stability which will allow you to take larger time steps without getting unrealistic results (although beware that “realistic” does not imply “accurate”).

Implicit finite difference stencil

The idea is the same, but our “T-shaped stencil” is now flipped upside down, and we get, 1k(un+1munm)=1h2(un+1m12un+1m+un+1m+1)unm=kh2un+1m1+(1+2kh2)un+1mkh2un+1m+1.

We can write it as a product of matrices to get un=Bun+1, where B has the same form as B in (???) but with a=1+2kh2 and b=kh2

Then, at each step we need to solve Equation (???) for un so: un+1=B1un

Adding in the reaction term

Let us now see how to update our value of u, when we are looking at a Reaction-Diffusion system like, ut=D2ux2+f(u).

Recall that if we just had, ut=f(u), then we simulate it using the forward Euler method so, un+1=f(un)k. where k=Δt.

Therefore, for the explicit method, we combine (???) and (???) to obtain that the update at each time step is, un+1=Bun+f(un)k

If we use the implicit method for the diffusion, and the explicit method, then we need to solve Bun+1=un+f(un)k, for un+1. As a result we get that at each step we need to compute, un+1=B1(un+f(un)k).

Bonus: Matlab Code

First, take a look at the following screencast. You can find all of the relevant Matlab code below so that you can play around with it on your own.

This code was used in the screencast above.

  1. %%%%%%%%%%%%%%%%%%%%%%%
  2. %%% Screencast code %%%
  3. %%%%%%%%%%%%%%%%%%%%%%%
  4.  
  5.  
  6. %%% Forward Finite difference method %%%
  7.  
  8. % At each time step compute: u(n+1) = B u(n)
  9. % How do we create B in MATLAB?
  10.  
  11. % Parameters
  12. D = 2; % Diffusion constant
  13. dx = 0.2; dt = 0.1; % space and time step
  14. Nx = 6; % Number of grid points
  15.  
  16. % Making the Matrix
  17. a = (1-2*D*dt/dx^2); % Diagonal values
  18. b = D*dt/dx^2; % Off-diagonal values
  19. main = a*sparse(ones(Nx,1));
  20. off = b*sparse(ones(Nx-1,1));
  21. B = diag(main) + diag(off,1) + diag(off,-1); % Make B
  22. % Apply periodic boundary conditions
  23. B(1, end-1) = b;
  24. B(end, 2) = b;
  25.  
  26.  
  27. % Useful MATLAB functions we used
  28. % diag
  29. % imagesc
  30. % sparse

This code will create the plots that were used to make the animations shown in Turing II.

  1. % The Mathematics of Patterns 2013
  2. % Turing Instabilities II: Gierer Meinhardt system in 1D
  3.  
  4. % This MATLAB code will show you how to numerically simulate the 1D Gierer
  5. % Meinhardt system, and make the (simplified) version of the animations in
  6. % the Turing Instabilities II section of the website.
  7. clear all
  8. close all
  9.  
  10. %%%%%%%%%%%%%%%%%%%%%
  11. %%% Set-up %%%
  12. %%%%%%%%%%%%%%%%%%%%%
  13.  
  14. % Parameter values
  15. bc = 0.35;
  16. Du = 1; Dv = 30; % Diffusion constants
  17.  
  18. % Grid and initial data:
  19. % w = 10; %no pattern
  20. w = 80; % pattern
  21.  
  22. Nx = 500; % How many points we want to discretize our domain with
  23. x = linspace(0,w, Nx);
  24. dx = x(2) - x(1);
  25.  
  26. dt = 1; % size of our time step
  27. t = 0:dt:400;
  28. Nt = length(t); % Number of time points
  29.  
  30. % Set up for the surface
  31. [X, T] = meshgrid(x, t);
  32. U = 0*X;
  33. V = 0*X;
  34.  
  35. % Easier to deal with column vectors
  36. x = x(:);
  37. t = t(:);
  38.  
  39. %Initial conditions: small perturbation away from steady state
  40. u = 1/bc*ones(length(x),1) + 0.01*rand(Nx, 1);
  41. v = 1/bc^2*ones(length(x),1);
  42.  
  43. % Save initial conditions
  44. U(1,:) = u;
  45. V(1,:) = v;
  46.  
  47.  
  48. %%%%%%%%%%%%%%%%%%%%%%%%%
  49. %%% Making the matrix %%%
  50. %%%%%%%%%%%%%%%%%%%%%%%%%
  51.  
  52. % To begin, let us recall how to set up the matrices used in the explicit
  53. % and implicit finite difference methods.
  54.  
  55.  
  56. %%% Forward (explicit) method %%%
  57. % We want a tridiagonal matrix (see notes for details)
  58. a = (1-2*Du*dt/dx^2); % values along the diagonal
  59. b = Du*dt/dx^2; % values in the off-diagonal
  60. main = a*sparse(ones(Nx,1));
  61. off = b*sparse(ones(Nx-1,1));
  62. Bu = diag(main) + diag(off,1) + diag(off,-1); %Still a sparse matrix
  63. % Satisfying boundary conditions
  64. Bu(1, end-1) = b;
  65. Bu(end, 2) = b;
  66. % To have a more numerically stable code, we use the implicit method.
  67.  
  68. %%% Backward (implicit) method %%%
  69. % For u
  70. a = (1+2*Du*dt/dx^2); % values along the diagonal
  71. b = Du*dt/dx^2; % values in the off-diagonal
  72. main = a*sparse(ones(Nx,1));
  73. off = -b*sparse(ones(Nx-1,1));
  74. Bu = diag(main) + diag(off,1) + diag(off,-1); %Still a sparse matrix
  75. % Satisfying boundary conditions
  76. Bu(1, end-1) = -b;
  77. Bu(end, 2) = -b;
  78.  
  79. % Same thing for v
  80. a = (1+2*Dv*dt/dx^2); b = Dv*dt/dx^2;
  81. main = a*sparse(ones(Nx,1));
  82. off = -b*sparse(ones(Nx-1,1));
  83. Bv = diag(main) + diag(off,1) + diag(off,-1);
  84. Bv(1, end-1) = -b;
  85. Bv(end, 2) = -b;
  86.  
  87. %%%%%%%%%%%%%%%%%%%%%%%%%
  88. %%% Plotting %%%
  89. %%%%%%%%%%%%%%%%%%%%%%%%%
  90.  
  91. figure(1); %create new figure
  92. plot(x,u,'g.-', 'linewidth',1);
  93. hold on;
  94. plot(x,v,'r.-', 'linewidth',1);
  95. hold off;
  96.  
  97. axis([-1 80 -.01 15.01]) % Fix axis limits
  98.  
  99. for j = 1:Nt
  100. % f and g are the reaction terms in the G-M system
  101. f = u.^2./v-bc*u;
  102. g = u.^2 - v;
  103.  
  104. % At each step we need to solve the system
  105. u = Bu\(u + dt*f); % backward Euler
  106. v = Bv\(v + dt*g);
  107. % Plot
  108. plot(x,u,'g.-', 'linewidth',1);
  109. hold on;
  110. plot(x,v,'r.-', 'linewidth',1);
  111. hold off;
  112. axis([-1 80 -.01 15.01])
  113. title(['t = ', num2str(j*dt)],'fontsize',24)
  114. drawnow;
  115. % Save for surface
  116. U(j,:) = u;
  117. V(j,:) = v;
  118. end
  119.  
  120.  
  121. %%%% Plotting the surface %%%%
  122. figure(2);
  123. s = surf(x, t, U)
  124. set(s, 'EdgeColor', 'none', 'FaceColor', 'interp'); % Sets up the colors
  125. xlabel('x')
  126. ylabel('t')
  127. zlabel('u')
  128.  
  129. %%%% contour plot %%%
  130. figure(3);
  131. p = pcolor(x, t, U);
  132. set(p, 'EdgeColor', 'none', 'FaceColor', 'interp');