抛物型方程的差分方程

数值解

---

例1：用Crank-Nicolson格式和追赶法求解抛物型方程

求解一维抛物方程的初边值问题：

$$\begin{eqnarray*} && \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2}, \ 0<x<\pi, t>0,\\ && u(0,t)=u(\pi,t)=0, \ t>0,\\ && u(x,0)=sin(x), \ 0<x<\pi \end{eqnarray*}$$
解：

clc
clear
xb = [0, pi]; %boundary of x
m = input('input the grid number of x:');
n = input('input the grid number of t:');
h = (xb(2)-xb(1))/m; %step value of x
k = 1/n; %step value of time
r = k/h^2; %ratio of grid
u = zeros(n+1, m+1); % the value of every grid
x = h*((1:m+1)-1);
u(1, :) = sin(x); %initial
u(:, 1) = 0; %boundary condition
u(:, end) = 0;
%A = zeros(m-2, 3); %the coefficient matrix of equation
d = zeros(1, m-1);
g = zeros(1, m-1);
w = zeros(1, m-1);
alpha = r/2;
beta = 1+r;
gamma = r/2;
for N = 2:n+1
    for M = 2:m
        ite = M-1;
        d(ite) = r/2*(u(N-1, M-1)+u(N-1, M+1))+(1-r)*u(N-1, M);
        if ite==1 %zhui
            g(1) = d(1)/beta;
            w(1) = gamma/beta;
        else
            g(ite) = (d(ite)+alpha*g(ite-1))/(beta-alpha*w(ite-1));
            w(ite) = gamma/(beta-alpha*w(ite-1));
        end
    end
    for M = m:-1:2 %gan
        ite = M-1;
        if M==m
            u(N, M) = g(ite);
        else
            u(N, M) = g(ite)+w(ite)*u(N, M+1);
        end
    end
end
X = xb(1):h:xb(2);
T = 0:k:n*k;
[X, T] = meshgrid(X, T);
surf(X, T, u);