向前差分算例

利用差分法求解下列问题：

??u?2u?2?2,0?x?1,t?0,??t?x??u(0,t)?u(1,t)?0,t?0??u(x,0)?sin(?x)?x(1?x),0?x?1???

x方向dx?0.1，t方向上dt?0.01。

在t?0.25观察数值解与精确解u?e??tsin(?x)?x(1?x)的误差。

一、算法描述 1. 网格剖分

取x??0,1?,t??0,0.25?，

2x方向dx?0.1，t方向上dt?0.01。

2. 差分格式

j?1j??u?ui?ui， ???dt??t?ij??2u?uij?1?uij?uij?1， ?2??2?xdx??ijr?dtdx2，

uij?1?r*uij?1?(1?2r)uij?r*uij?1?2*dt。

3. 初边值处理

可以先设置21?100的0矩阵，然后可以直接给出初值条件，边界条件为0可以不用设置。

二、程序：

%***********清理内存************************************ clear all clc

%**********网格划分************************************** N=10;%总空间长

T=25;%求解时刻或总的时间步 dx=0.1;%空间步长 dt=0.01;%时间步长

%***********定义0矩阵************************************ u=zeros(N+1,T+1);

%************边界条件************************************ u(1,1:T+1)=0;

u(N+1,1:T+1)=0;

%************初始条件************************************ for i=1:N+1 x=(i-1)*dx;

u(i,1)=sin(pi*x)+x.*(1-x);%得到空间上第1层的值 end

%************差分迭代************************************* r=dt/dx^2; for j=1:T

for i=2:N

u(i,j+1)=r*u(i+1,j)+(1-2*r)*u(i,j)+r*u(i-1,j)+2*dt; %由第1层的值得到2层及其以上的空间上的值 end end

disp(u');%在屏幕上显示全部求解结果 U=u(1:N+1,T+1);%提取t=0.25时刻的值

%*********绘制0.25时刻的数值解与精确解图******************** x=zeros(1,N+1); for i=1:N+1

x(i)=(i-1)*dx; end

plot(x,U','o') hold

u_xt=exp(-pi*pi*T*dt).*sin(pi.*x)+x.*(1-x); plot(x,u_xt,'r')

e=u_xt-U' %求解误差

%**********************************************************

三、结果

注：红色曲线代表精确解，圆圈代表数值解。

本文来源：https://www.bwwdw.com/article/g7dh.html