偏微分方程数值解课程设计

课程设计报告
课程：偏微分方程数值解学号：
姓名：
班级：
教师：
《偏微分方程数值解》
课程设计指导书
一.课程设计的目的
1.帮助掌握偏微分方程数值解相关知识。

2.理解偏微分方程数值解差分隐格式解决自由振动方程问题的方法。

3.锻炼编写程序代码的能力。

二.设计名称
差分法求自由振动问题的周期解。

三.设计要求
1.要求写出差分隐格式的理论方法。

2.要求编写matlab 程序，画出函数图形。

3.要求写出实验总结及心得体会。

四.设计题目
用差分法求自由振动问题的周期解：
2222000,,0|0,|sin (0,)(2,)t t u u
x t t x u u x t u t u t π==⎧∂∂-=-∞<<∞>⎪∂∂⎪
∂⎪==⎨∂⎪
=⎪⎪⎩
要求用差分隐格式求解，其中14
θ=。

五.设计细则 1.区域剖分：
构造上式的差分逼近，取空间步长h 和时间步长τ，用两族平行直线
⎩⎨
⎧===±±===Λ
Λ,2,1,0,,
,2,1,0,n n t t j jh x x n j τ 作矩形网格。

2.离散格式：
显格式：
于网点),(n j t x 用Taylor 展式，并整理方程得：
⎪⎪⎩
⎪⎪
⎨⎧--++=+-++==-+-++-121121102
10102100
)1(2)(),()()1()]()([2),(n j n j n j n j n j j j j j j j j u u r u u r u x x r x x r u x u τϕϕϕϕϕ
隐格式：
上述显格式并不是绝对稳定的差分格式，为了得到绝对稳定的差分格式，用第1-n 层、
n 层、1+n 层的中心差商的权平均去逼近xx u ，得到下列差分格式：
⎪⎪⎪⎩
⎪⎪
⎪⎨⎧+-++--++-=+-+-++==----+-++-+++-++-]22)21(2[2),
()()1()]()([2),(2
111112112111112
211102
10102100h u u u h u u u h u u u a u u u x x r x x r u x u n j n j n j n j n j n j n j n j n j n j n j n j j j j j j j j θθθττϕϕϕϕϕ其中10≤≤θ是参数。

当0=θ时就是显格式，而当4
1
=θ时可以证明该格式绝对稳定。

隐格式的矩阵形式是：
⎥⎥⎥
⎥⎥⎥⎥⎥⎥⎦⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡-+-+-+-+--+-+-+++122111121121
12222
222
2222221212121J J j n J n J n j n n z z z z z u u u u u r r r r r r r r r r r r M M M M M M M θθθθ
θθθθθ
θ
θθ
其中：
1
111111122]2()2)(21[(-----+-+-++-++--=n j n j n j n j n j n j n j n j j u u u u u u u u r z θθ
3.格式稳定性： 1）显格式： 显格式稳定的充分必要条件是：网格比1<r 。

2）隐格式：
当4
1
=
θ时隐格式绝对稳定。

4.数值例子：
可以证明 )
(),(t x e
t x u +=是波动方程222
2u u
t x
∂∂=∂∂ 的一个解析解。

那么，为了更精确的得到误差估计，在这里选取)
(),(t x e t x u +=作数值实验。

取0201x t π
≤≤≤≤，并且将时间步长10等分，空间步长100等分（即
1/10,2/100h τπ==）。

这样网格比/15
r ah π
τ==<，从稳定性分析可知，此时格式稳
定。

六.程序代码： 1.主函数：
T=1;%取时间长度为1 b=0.5;
h=2*pi/100; tao=1/10;
f=inline('0','x','t');%f=0 fx1=inline('0');
fx2=inline('sin(x)');
ft1=inline('sin(t)');%此题中取ft1=ft2=sin(t)=sin(2*pi+t) ft2=inline('sin(2*pi+t)');
[X,Y,U,Z]=chfenmethed(f,fx1,fx2,ft1,ft2,T,h,tao); mesh(X,Y,U);%绘制函数图像 shading flat;
xlabel('X','FontSize',14); ylabel('t','FontSize',14); zlabel('U','FontSize',14); title('函数'); figure;
mesh(X,Y,Z);%绘制误差图像 xlabel('X','FontSize',14); ylabel('t','FontSize',14); zlabel('Z','FontSize',14); title('误差');
2.编写差分隐格式函数：
function [X,T,U,Z]=chfenmethed(f,fx1,fx2,ft1,ft2,T,h,tao) %u_tt-u_xx=f(x,t) 0<x<2*pi,0<t<T
%u(0,t)=ft1,u(2*pi,t)=ft2,此题中ft1=ft2 %u(x,0)=fx1(x),此题中fx1=0
%u_t(x,0)=fx2(x),此题中fx2=sin(x) x=0:h:2*pi; t=0:tao:T; m=length(x); n=length(t); s=tao/h;
[X,T]=meshgrid(x,t); Z=zeros(n,m); U=zeros(n,m);
for i=2:m-1
U(1,i)=feval(fx1,x(i));
U(2,i)=U(1,i)+tao*feval(fx2,x(i))+tao^2/2*(1/h^2* ...
(feval(fx1,x(i+1))-2*feval(fx1,x(i))+feval(fx1,x(i-1))+feval(f,x(i),0))); Z(2,i)=abs(U(2,i)-f0(x(i),t(2)));
end
for j=1:n
U(j,1)=feval(ft1,t(j));
U(j,m)=feval(ft2,t(j));
end
A=-0.5*s^2*ones(1,m-2);
C=A;
B=(1+s^2)*ones(1,m-2);
UU=zeros(1,m-2);
f1=UU;
for i=3:n
for j=2:m-1
UU(j-1)=f0(x(j),t(i));
f1(j-1)=0.5*s^2*U(i-2,j-1)-(1+s^2)*U(i-2,j)...
+0.5*s^2*U(i-2,j+1)+2*U(i-1,j)...
+tao^2*feval(f,x(j),t(i-1));
end
f1(1)=f1(1)+0.5*s^2*U(i,1);
f1(end)=f1(end)+0.5*s^2*U(i,m);
U(i,2:m-1)=zgf(A,B,C,f1);
Z(i,2:m-1)=abs(U(i,2:m-1)-UU);
end
3.编写迭代函数，用以实现差分隐格式：
function x=zgf(A,B,C,f)
n=length(B);
B1=zeros(1,n-1);
Y=zeros(1,n);
x1=zeros(1,n);
B1(1)=C(1)/B(1);
for i=2:n-1
B1(i)=C(i)/(B(i)-A(i)*B1(i-1));
end
Y(1)=f(1)/B(1);
for i=2:n
Y(i)=(f(i)-A(i)*Y(i-1))/(B(i)-A(i)*B1(i-1));
end
x1(n)=Y(n);
for i=n-1:-1:1
x1(i)=Y(i)-B1(i)*x1(i+1); end
x=x1;
4.编写精确函数用以求误差：
function z=f0(x,t)%精确解函数
z=exp(x+t);
七.设计结果与分析：
最终得到的两张图像如下：
八.设计体会与建议：
设计成绩：教师签名：
年月日。