Scientific Computing

Scientific Computing Partial Differential Equations Explicit Solution of Heat Equation

The one-dimensional Heat Equation for the temperature u(x,t) in a rod at time t is given by: where c >0, T>0 are constants One Dimensional Heat Equation

We will solve this equation for x and t values on a grid in x-t space: One Dimensional Heat Equation Approximate Solution uij=u(xi, tj) at grid points

To approximate the solution we use the finite difference approximation for the two derivatives in the heat equation. We use the forward-difference formula for the time derivative and the centered-difference formula for the space derivative: Finite Difference Approximation

Then the heat equation (ut =cuxx ) can be approximated as Or, Let r = (ck/h2) Solving for ui,j+1 we get: Finite Difference Approximation

Note how this formula uses info from the j-th time step to approximate the (j+1)-st time step: One Dimensional Heat Equation

On the left edge of this grid, u0,j = g0,j = g0(tj). On the right edge of this grid, un,j = g1,j = g1(tj). On the bottom row of the grid, ui,0 = fi = f(xi). Thus, the algorithm for finding the (explicit) solution to the heat equation is: One Dimensional Heat Equation

function z = simpleHeat(f, g0, g1, T, n, m, c) %Simple Explicit solution of heat equation h = 1/n; k = T/m; r = c*k/h^2; % Constants x = 0:h:1; t = 0:k:T; % x and t vectors % Boundary conditions u(1:n+1, 1) = f(x)'; % Transpose, since it’s a row vector u(1, 1:m+1) = g0(t); u(n+1, 1:m+1) = g1(t); % compute solution forward in time for j = 1:m u(2:n,j+1) = r*u(1:n-1,j) + (1-2*r)*u(2:n,j) + r*u(3:n+1,j); end z=u'; mesh(x,t,z); % plot solution in 3-d end Matlab Implementation

Usage: f = inline(‘x.^4’); g0 = inline(‘0*t’); g1 = inline(‘t.^0’); n=5; m=5; c=1; T=0.1; z = simpleHeat(f, g0, g1, T, n, m, c); Matlab Implementation

Calculated solution appears correct: Matlab Implementation

Try different T value: T=0.5 Values seem chaotic Matlab Implementation

Why this chaotic behavior? Matlab Implementation

Matrix Form of Solution

This is of the form Start: u(0)=u(:,0)=f(:) Iterate: Matrix Form of Solution

Now, suppose that we had an error in the initial value for u(0), say the initial u value was u(0)+e, where e is a small error. Then, under the iteration, the error will grow like Ame. For the error to stay bounded, we must have Thus, the largest eigenvalue of A must be <= 1. Matrix Form of Solution

Let A be a square nxn matrix. Around every element aii on the diagonal of the matrix draw a circle with radius Such circles are known as Gerschgorin disks. Theorem: Every eigenvalue of A lies in one of these Gerschgorin disks. Gerschgorin Theorem

Example: The circles that bound the Eigenvalues are: C1: Center point (4,0) with radius r1 = |2|+|3|=5 C2: Center point (-5,0) with radius r2=|-2|+|8|=10 C3: Center Point (3,0) with radius r3=|1|+|0|=1 Gerschgorin Theorem

Example: Actual eigenvalues in red Gerschgorin Theorem

Example: C1: Center point (1,0) radius r1 = |0|+|7|=7 C2: Center point (-5,0) radius r2=|2|+|0|=2 C3: Center Point (-3,0) radius r3=|4|+|4|=8 Gerschgorin Theorem

Circles: center=(1-2r), radius = r or 2r. Take largest radius 2r. Then, if λ is an eigenvalue of A, we have Matrix Form of Solution 1-2r 0

So, the eigenvalue satisfies: For the error to stay bounded, we need Thus, we need For our first example, we had T=0.1 and so, we expect the solution to be stable. For the second example, T=0.5, we have r =2.5, so the error in the iterates will grow in an unbounded fashion, which we could see in the numerical solution. Matrix Form of Solution

Scientific Computing