@why-math
2017-05-17T11:06:57.000000Z
字数 1893
阅读 4134
teaching
参考教材: 微分方程数值解(第四版),李荣华,刘播;
给定如下求解个维热传导方程的数值格式:
运行如下测试脚本 errortest.m
:
%% 针对不同的 r , 测试向前差分格式中初始误差随迭代的传播情况
%
% $$ u_j^{n+1}=ru_{j+1}^n+(1-2r)u_j^n+ru_{j-1}^n+\tau f_j $$
%
% 作者:魏华祎 <weihuayi@xtu.edu.cn>
N = 10;
epsilon = 0.0001;
maxit = 100;
r = 5/9;
E1 = errorvary(N, r, epsilon, maxit);
r = 5/11;
E2 = errorvary(N,r,epsilon, maxit);
% 显示 r= 5/9 初始误差和第 100 个时间步的误差
disp(E1(:,[1,100]));
% 显示 r= 5/11 初始误差和第 100 个时间步的误差
disp(E2(:,[1,100]));
函数 errorvary
定义如下:
function E = errorvary(N,r,epsilon,maxit)
%% ERRORVARY 给定有 N 个节点的网格,实验检验不同 r 清况下,初始误差随时间推进的变化情况
%
% 输入参数:
% N 网格节点数
% r 网格比
% epsilon 第 round(N/2) 个网格节点处的初始误差
% maxit 最大迭代步数
% 输出参数:
% E N*maxit 阶矩阵,E(:,i) 表示第 i 个时间层上的每个网格点处的误差
%
% 作者:魏华祎 <weihuayi@xtu.edu.cn>
E = zeros(N,maxit);
E(round(N/2),1) = epsilon;
d = 1 - 2*ones(N-2,1)*r;
c = ones(N-3,1)*r;
A = diag(c,-1) + diag(c,1)+diag(d);
for i = 2:maxit
E(2:end-1,i) = A*E(2:end-1,i-1);
end
end
其中
空间区间为;
时间区间;
;
真解
把真解代入上面模型, 即可得到相关的参数.
最大模误差定义如下:
即所有网格点处数值解和真解误差绝对值的最大值, 最大模误差 与时间步长 和空间步长 满足如下关系:
所以,当 变为 , 变为 时, 应变为 .
示例Matlab代码请参考这里;
用 Ritz-Galerkin 方法求边值问题