[微]中心差分格式求解两点边值问题 参考程序

本问题的难度比较低，主要是求解一个三对角线性方程组。该线性方程组的求解使用了函数TridiagAXb。下面的程序是主程序，后面会给出函数TridiagAXb的代码。

代码:


% 本程序使用中心差分格式计算两点边值问题
% -u'' = f, u(0) = u(1) = 0
% u_exact = 1 - (1-exp(-10))*x - exp(-10*x)
% f = 100*exp(-10*x)
%  n        误差结果-无穷范数 err_inf
%  20        1.3742903775545e-2
%  40        3.476685596693e-3
%  80        8.712085690456e-4
%  160      2.179740562154e-4
% 2阶收敛速度
%
%  n        误差结果-无穷范数 err_L2
%  20        9.277264457428e-3
%  40        2.341306410922e-3
%  80        5.867024245239e-4
%  160      1.467616112220e-4
% 2阶收敛速度


clear; close all;format long;
n = 160;
h = 1/n;
x = h:h:1-h;
u_exact = zeros(n-1,1);
u = zeros(n-1,1);
b = zeros(n-1,1);
for i = 1:n-1
    b(i) = f(x(i))*h^2;
    u_exact(i) = 1 - (1-exp(-10))*x(i) - exp(-10*x(i));
end
u = TridiagAXb(-1,2,-1,b);

% 画图
plot(x,u,'+',x,u_exact,'r')
% 误差计算
err_inf = max(abs(u - u_exact))
err_L2 = 0;
for i = 1:n-1
    err_L2 = err_L2 + h*(u(i) - u_exact(i))^2;
end
err_L2 = sqrt(err_L2)

函数TridiagAXb

代码:


% 三对角方程求解函数  M*M
% AX = B,  A的主对角线元素为b，右上次对角线元素为a，左下次对角线元素为c
function X = TridiagAXb(a,b,c,B)
M = length(B);
X = zeros(M,1);
D = zeros(1,M);
D(1) = b;
for k = 1:M-1
    temp = c/D(k);
    D(k+1) = b - a*temp;
    B(k+1) = B(k+1) - B(k)*temp;
end
X(M) = B(M)/D(M);
for k = M-1:-1:1
    X(k) = (B(k) - X(k+1)*a)/D(k);
end
return;