数值微分方程第二次上机实验Word格式文档下载.docx
- 文档编号:15714095
- 上传时间:2022-11-15
- 格式:DOCX
- 页数:14
- 大小:335.84KB
数值微分方程第二次上机实验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.由于系数矩阵的特殊性质,可以使用共轭梯度发进行求解,在以后会进行尝试。
【结论】
通过实验结果可知:
五点差分格式法求解椭圆型方程的精度较高,耗时较少,得到的结果比较好;
而利用数值方法求解,当维数变大后,虽然结果精度依然较高,但所耗时间明显增加.
指导教师评语及成绩:
成绩:
指导教师签名:
批阅日期:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 微分方程 第二次 上机 实验