数值微分方程第二次上机实验Word格式文档下载.docx
《数值微分方程第二次上机实验Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数值微分方程第二次上机实验Word格式文档下载.docx（14页珍藏版）》请在冰豆网上搜索。
Windows8专业版
实验内容：
【实验方案设计】
第一步，将书上关于三种插值方法的内容转化成程序语言，用MATLAB实现；
第二步，分别运用以上求解初值问题的程序求解具体问题。
【实验过程】
（实验步骤、记录、数据、分析）
主要步骤：
首先，分析问题，根据分析设计MATLAB程序；
其次，利用程序算出问题答案；
最后，对答案进行分析，并得出结论。
实验:
用五点差分格式求解课本第101页第4题的差分方程组，并用对称矩阵与向量表示。
要求如下:
（1）边值条件尝试：
常数（10）或者函数（10*sin（x））。
（2）步长尝试：
h=（0.5）^n,n=2,3,5。
（3）Au=g，用内置方法和任一数值解法。
（1）边界条件为常数时,编写如下程序:
function[A,g,U]=pde（h）
x=0:
h:
1;
n=length（x）;
nx=length（x）-1;
U=zeros（n）;
U（n,:
）=0;
A=zeros（（nx）*n）;
g=zeros（nx*n,1）;
B=eye（nx）*4;
fori=1:
nx-1
B（i+1,i）=-0.5;
B（i,i+1）=-0.5;
end
B（1,1）=2;
K=eye（nx）*（-3）;
K（1,1）=-1.5;
A（1:
nx,1:
nx）=B;
nx,nx+1:
2*nx）=K;
gg=ones（nx,1）*8*h*h;
gg
（1）=4*h*h;
g（1:
nx）=gg;
n-2
A（i*nx+1:
（i+1）*nx,（i-1）*nx+1:
i*nx）=K;
（i+1）*nx,（i+1）*nx+1:
（i+2）*nx）=K;
（i+1）*nx,i*nx+1:
（i+1）*nx）=2*B;
g（i*nx+1:
（i+1）*nx）=gg;
A（（n-1）*nx+1:
n*nx,（n-2）*nx+1:
（n-1）*nx）=K;
C=eye（nx）*4;
C（i+1,i）=-0.5;
C（i,i+1）=-0.5;
C（1,1）=2;
n*nx,（n-1）*nx+1:
n*nx）=C;
gh=ones（nx,1）*（8*h*h-3*h*10）;
gh
（1）=4*h*h-1.5*h*10;
g（（n-1）*nx+1:
n*nx）=gh;
u=A\g;
n
U（1:
n-1,i）=u（（i-1）*nx+1:
i*nx）;
mesh（U）;
运行程序：
取h=1/4,在command命令框中输入如下命令:
>
[A,g,U]=pde（1/4）
得到如下结果:
A矩阵为：
g=
0.2500
0.5000
-3.5000
-7.0000
U=
-8.2584-8.6250-9.5628-11.0834-13.1994
-7.6585-8.0182-8.9387-10.4402-12.5473
-5.9000-6.2294-7.0712-8.4804-10.5380
-3.1651-3.4035-4.0011-5.0757-6.8741
00000
分别取h=（1/2）^n,n=2,3,5,只列出所做图像:
（2）上边界条件取函数，10*sin（x）
则包含上边界条件部分程序改为
fori=n/h:
n/h:
1
gh=ones（nx,1）*（8*h*h-3*h*10*sin（i））;
gh
（1）=4*h*h-1.5*h*10*sin（0）;
运行结果如下：
4.70444.63584.59784.59484.6390
4.41024.33994.29524.27404.2716
3.53743.46763.42193.39703.3895
2.08402.02241.98661.96791.9621
分别取h=（1/2）^n,n=2,3,5,只列出所做图像
（3）取上边界条件为-u，则包含上边界点的相关程序改为
C=eye（nx）*（4-3*h）;
C（1,1）=2+1.5*h;
n*nx）=gg;
（4）使用数值解法求解
i.使用LU列主元分解，后求解
编写程序如下:
function[L,U,u]=Lu2（A）
n=length（A）;
L=zeros（n）;
B=zeros（1,n）;
s=0;
n-1
t=1;
[s,t]=max（abs（A（i:
n,i）））;
m=t+i-1;
u（i）=m;
B=A（m,:
）;
A（m,:
）=A（i,:
A（i,:
）=B;
ifA（i,i）~=0
A（i+1:
n,i）=A（i+1:
n,i）/A（i,i）;
n,i+1:
n）=A（i+1:
n）-A（i+1:
n,i）*A（i,i+1:
n）;
else
stop;
end
L=tril（A）;
L（1:
n+1:
end）=1;
U=（tril（A'
））'
;
u（n）=n;
function[x]=Solve（A,b）
[L,U,u]=Lu2_lx（A）;
temp=b（i）;
b（i）=b（u（i））;
b（u（i））=temp;
x=（L*U）\b;
使用该方法对h=1/4进行运算，比较时间与准确度:
tic,[x]=Solve（A,g）;
toc
时间已过0.003258秒。
norm（x-u）
ans=
1.0582e-11
使用该方法对h=1/32进行运算，比较时间与准确度
toc
时间已过4.695820秒。
norm（x-q）
5.3498e-13
可得该方法准确度还不错，但是维数大时，耗时较长。
ii.使用sor迭代
function[Y,k]=Sor（A,b,w,a）
k=0;
x=ones（1,n）;
x=x'
while（norm（A*x-b,2）>
=a）
fori=1:
x（i）=w*（b（i）+A（i,1:
i-1）*x（1:
i-1）+A（i,i+1:
n）*x（i+1:
n））/A（i,i）+（1-w）*x（i）;
k=k+1;
Y=x;
End
在command命令框中输入如下命令:
tic,[Y,k]=Sor（A,g,w,a）;
时间已过54.506909秒。
iii.由于系数矩阵的特殊性质，可以使用共轭梯度发进行求解，在以后会进行尝试。
【结论】
通过实验结果可知:
五点差分格式法求解椭圆型方程的精度较高,耗时较少,得到的结果比较好;
而利用数值方法求解,当维数变大后,虽然结果精度依然较高,但所耗时间明显增加.
指导教师评语及成绩：
成绩：
指导教师签名：
批阅日期：
升级会员 