偏微分方程数值实验报告(算法及matlab实现)(2) 下载本文

偏微分方程实验报告二

实验题目:利用向前差分求解如下一维热传导方程,并观察最大规模误差变化情况:

?u?2u?a2?f(x,t),?t?xu(0,t)?uL(x),u(1,t)?uR(x), u(x,0)?u0(x),x?[0,1],t?[0,0.1].其中a?1;真解为u(x,t)?sin(2?x)e10t 最大模误差定义如下:

E?maxu(xi,tj)?U(i,j)

xi,tj

实现算法:对于初边值问题(抛物方程)

?u?2u?a2?f(x),(x,t)?G,?t?x (1) u(x,0)??(x),0?x?l,u(0,t)?u(l,t)?0,0?t?T,其相应的有限差分法的算法如下:

1.对求解区域做网格剖分,得到计算网格.取空间步长和时间步长为h?1T,??,分别对NM空间变量X所属的区间[0,l]和时间变量t所属的区间[0,T]做如下均匀剖分:

0?x0?x1?...?xN?l,0?t0?t1?...tM?T

其中xj?jh,tk?k?.

用两族平行直线x?xj(0,1,...,N)和t?tk(k?0,1,...,M)将矩形域G分割成矩形网格. 2.对微分方程中的各阶导数进行差分离散,得到差分方程.用uj表示差分解在网格节点处的分量,下面采用逐层计算的思想建立差分格式. (xj,tk)设从第0个时间层到第k?0个时间层的差分解分量

kuij,i?0,...,k;j?0,...,N

已经求得.

下面建立第k?1个时间层(简称当前层)上的向前差分方程,对当前层上的任意内节点

(xj,tk?1),规定其差分格式涉及的模板点为

对处的微分方程 (xj,tk)?u?2u?a2?t(xj,tk)?x(x中的导数用差商近似:

?f(xj) (1)

j,tk)?u(xj,tk)?t?2u(xj,tk)?x2?u(xj,tk?1)?u(xj,tk)(向前差商)??u(xj?1,tk)?2u(xj,tk)?u(xj?1,tk)h2(2)

(二阶中心差商)记fj?f(xj),将(2)代入(1),则有关于的向前差分格式 (xj,tk?1)?1uk?ukjj?在(3)两边同乘?,整理可得

?akkukj?1?2uj?uj?1h2?fj (3)

?1kkuk?rukjj?1?(1?2r)uj?ruj?1??fj (4)

其中,r?a?h2称为网格比(简称网比).

3.根据初边值条件,进行初边值处理.由(4)可得

?1kkuk?ruk;k?0,...,M?1jj?1?(1?2r)uj?ruj?1??fj,j?1,...,N?1u0j?u(xj,0)??(xj),j?1,...,N?1kku0?u(0,tk)?0,uN?u(l,tk)?0,k?0,1,...,M(5)

令N?1维(第k个网格层)差分解向量和右端向量

kkTTUk?(u1k,u2,...,uN?1),F?(f1,f2,...,fN?1)

则向前差分格式(4)的矩阵表示为:

Uk?1?A0Uk??F

其中

4.最后求解线性代数方程组,得到数值解向量Un?1

.

程序代码:

function [X,T,U]=heat_equation_fd1d(NS,NT,pde) %% HEAT_EQUATION_FD1D 利用有限差分方法计算一维热传导方程 % 输入参数:

% NS _x0012_整型,空间剖分段数 % NT _x0012_整型,时间剖分段数

% pde 结构体,待求解的微分方程模型的已知数据,如边界、初始、系数和右端项等条件 % 输出参数

% X 长度为 NS +1 的列向量,空间网格剖分 % T 长度为 NT +1 的行向量,时间网格剖分

% U (NS +1) *( NT +1) 矩阵,U(:,i)表示第i个时间层网格部分上的数值解 [X,h]=pde.space_grid(NS); [T,tau]=pde.time_grid(NT); N=length(X); M=length(T); r=pde.a()*tau/h/h; if r>=0.5

display('时间空间离散不满足向前差分的稳定条件'); end

U=zeros(N,M);

U(:,1)=pde.u_initial(X); U(1,:)=pde.u_left(T); U(end,:)=pde.u_right(T); %% 向前差分

d=1-2*ones(1,N-2)*r; c=ones(1,N-3)*r;

A=diag(c,-1)+diag(d)+diag(c,1); for i=2:M

RHS=tau*pde.f(X,T(i)); RHS(2)=RHS(2)+r*U(1,i-1);

RHS(end-1)=RHS(end-1)+r*U(end,i-1); U(2:end-1,i)=A*U(2:end-1,i-1)+RHS(2:end-1); end end

function emax=heat_equation_fd1d_error(X,T,U,u_exact)

ue=u_exact(X,T); %真解在网格节点处的值 emax=max(max(ue-U)); end

function showsolution (X,T,U) %% SHOWSOLUTION 以二元函数方式显示数值解 % 输入参数

% X 长度为 NS +1 的列向量,空间网格剖分N % T 长度为 NT +1 的行向量,时间网格剖分M

% U N*M 矩阵,U(:,i)表示第i个时间层网格部分上的数值解 [x,t] = meshgrid (X,T); mesh (x,t,U'); xlabel ('X'); ylabel ('T'); zlabel ('U(X,T)'); End

function showvarysolution (X,T,U) %% SHOWVARYSOLUTION 显示数值解随着时间的变化 % % 输入参数

% X 长度为 NS +1 的列向量,空间网格剖分N % T 长度为 NT +1 的行向量,时间网格剖分M

% U N*M 矩阵,U(:,i)表示第i个时间层网格部分上的数值解 M = size (U ,2) ; figure xlabel ('X'); ylabel ('U');

s = [X(1) ,X(end),min(min(U)),max(max(U))]; axis (s); for i = 1:M

plot (X,U(:,i)); axis (s); pause (0.0001) ;

title ([ 'T=',num2str(T(i)),'时刻的温度分布']) end

function u_exact=u_exact(X,T) u_exact=sin(2*pi*X)*exp(10*T); End

function pde=model_data()

pde=struct('u_initial',@u_initial,'u_left',@u_left,'u_right',@u_right,'f',@f,'time_grid',@time_grid,...

'space_grid',@space_grid,'a',@a); function [T,tau]=time_grid(NT) T=linspace(0,0.1,NT);%T为行向量 tau=0.1/(NT-1); end

function [X,h]=space_grid(NS) X=linspace(0,1,NS)';%X为列向量 h=1/(NS-1); end

function u=u_initial(X) u=sin(2*pi*X); end

function u=u_left(T)%左边值条件,x0位置不动,时间在变化 u=zeros(size(T)); end

function u=u_right(T)%右边值条件,xN位置不动,时间在变化 u=zeros(size(T)); end

function f=f(x,t)

f=10*exp(10*t)*sin(2*pi*x)+4*a*pi*pi*exp(10*t)*sin(2*pi*x); end

function a=a() a=1; end end

pde=model_data(); %模型数据结构体 main_test %向前差分

N=[11,21,31]; M=[101,201,301]; for i=1:3

[X,T,U]=heat_equation_fd1d(N(i),M(i),pde); showvarysolution(X,T,U); %以随时间变化方式显示数值解 showsolution(X,T,U);%以二元函数方式显示数值解

[emax]=heat_equation_fd1d_error(X,T,U,@u_exact); format short disp('emax'); disp([emax]); end

实验结果:在命令窗口输入>> main_test 回车可得运算结果和图像 >> main_test emax