# A Guide to MATLAB for Beginners and Experienced Users - Brian R.H.

**Download**(direct link)

**:**

**60**> 61 62 63 64 65 66 .. 91 >> Next

Generally speaking, it is best to understand some of the theory of partial differential equations before attempting a numerical solution like we have done here. However, for this particular case at least, the simple rule of thumb of keeping the coefficients of the iteration positive yields realistic results. A theoretical examination of the stability of this finite difference scheme for the one-dimensional heat equation shows that indeed any value of s between 0 and 0.5 will work, and it suggests that the best value of At to use for a given Ax is the one that makes s = 0.25. (See Partial Differential Equations: An Introduction, by Walter A. Strauss, John Wiley and Sons, 1992.) Notice that while we can get more accurate results in this case by reducing Ax, if we reduce it by a factor of 10 we must reduce At by a factor of 100 to compensate, making the computation take 1000 times as long and use 1000 times the memory!

Earlier we mentioned that the problem we solved numerically could also be solved analytically. The value of the numerical method is that it can be applied to similar partial differential equations for which an exact solution is not possible or at least not known. For example, consider the one-dimensional heat equation with a variable coefficient, representing an inhomogeneous material with varying thermal conductivity k(x),

For the first derivatives on the right-hand side, we use a symmetric finite difference approximation, so that our discrete approximation to the partial differential equations becomes

The Case of Variable Conductivity

At Ax" 2 Ax 2 Ax

where kj = k(a + j Ax). Then the time iteration for this method is

2A x 2A x

un+1 = sjj + un-1) + (1 - 2sj) un + 0.25 j - sj—1) j — u1]^),

n+1 —190 Chapter 9: Applications

where Sj = kjAt/(Ax)2. In the following M-file, which we called heatvc.m, we modify our previous M-file to incorporate this iteration.

function u = heatvc(k, x, t, init, bdry)

% Solve the 1D heat equation with variable coefficient k on % the rectangle described by vectors x and t with % u(x, t(1)) = init and Dirichlet boundary conditions % u(x(1), t) = bdry(1), u(x(end), t) = bdry(2).

J = length(x); N = length(t); dx = mean(diff(x)); dt = mean(diff(t)); s = k*dt/dx"2;

u = zeros(N,J); u(1,:) = init;

for n = 1:N-1

u(n+1, 2:J-1) = s(2:J-1).*(u(n, 3:J) + u(n, 1:J-2)) + ...

(1 - 2*s(2:J-1)).*u(n,2:J-1) + ...

0.25*(s(3:J) - s(1:J-2)).*(u(n, 3:J) - u(n, 1:J-2)); u(n+1, 1) = bdry(1); u(n+1, J) = bdry(2); end

Notice that k is now assumed to be a vector with the same length as x and that as a result so is s. This in turn requires that we use vectorized multiplication in the main iteration, which we have now split into three lines.

Let's use this M-file to solve the one-dimensional variable-coefficient heat equation with the same boundary and initial conditions as before, using k(x) = 1 + (x/5)2. Since the maximum value of k is 2, we can use the same values of At and Ax as before.

kvals = 1 + (xvals/5)."2;

uvals = heatvc(kvals, xvals, tvals, init, [15 25]); surf(xvals, tvals, uvals) xlabel x; ylabel t; zlabel uNumerical Solution of the Heat Equation

191

In this case the limiting temperature distribution is not linear; it has a steeper temperature gradient in the middle, where the thermal conductivity is lower. Again one could find the exact form of this limiting distribution, u(x, t) = 20(1 + (1/n)arctan(x/5)), by setting the t derivative to zero in the original equation and solving the resulting ordinary differential equation.

You can use the method of finite differences to solve the heat equation in two or three space dimensions as well. For this and other partial differential equations with time and two space dimensions, you can also use the PDE Toolbox, which implements the more sophisticated finite element method.

A SIMULINK Solution

We can also solve the heat equation using SIMULINK. To do this we continue to approximate the x derivatives with finite differences, but we think of the equation as a vector-valued ordinary differential equation, with t as the independent variable. SIMULINK solves the model using MATLAB's ODE solver, ode45. To illustrate how to do this, let's take the same example we started with, the case where k = 2 on the interval —5 < x < 5 from time 0 to time 4, using boundary temperatures 15 and 25, and initial temperature distribution of 15 for x < 0 and 25 for x > 0. We replace u(x, t) for fixed t by the vector u of values of u(x, t), with, say, x = -5:5. Here there are 11192 Chapter 9: Applications

values of X at which we are sampling u, but since u(x, t) is pre-determined at the endpoints, we can take u to be a 9-dimensional vector, and we just tack on the values at the endpoints when we're done. Since we're replacing ä2u/dX2 by its finite difference approximation and we've taken Ax = 1 for simplicity, our equation becomes the vector-valued ODE

ä u

~ç¿

= k( Au + c).

**60**> 61 62 63 64 65 66 .. 91 >> Next