数值分析 实验报告四.docx
- 文档编号:3140989
- 上传时间:2023-05-05
- 格式:DOCX
- 页数:8
- 大小:30.96KB
数值分析 实验报告四.docx
《数值分析 实验报告四.docx》由会员分享,可在线阅读,更多相关《数值分析 实验报告四.docx(8页珍藏版)》请在冰点文库上搜索。
数值分析实验报告四
数学与信息工程学院
实验报告
课程名称:
计算方法
实验室:
7404
实验台号:
班级:
11数学
(1)班
姓名:
陈露
实验日期:
2014年5 月21日
实验名称
线性方程组迭代解法及常微分方程数值解法
实验目的
和要求
1.Jacobi迭代法的Matlab实现;
2.G-S迭代法的Matlab实现;
3.选做:
超松弛迭代法。
4.Euler法的Matlab实现;
5.改进Euler法的Matlab实现
实验内容和步骤:
实验内容:
1、Thelinearsystem
hasthesolution
Usethejacobimethodwith
toapproximatethesolutiontothelinearsystemtowithin
inthe
norm.
2、ProgramtheG-Smethodandtestitontheseexamples:
.
3、给定初值问题
(1)用欧拉法(h=0.05)求数值解。
(2)用改进欧拉法(h=0.05)求数值解。
实验步骤:
Jacobi迭代法
1)实验编程
1.新建M文件JacobiIteration.m输入:
functionx=JacobiIteration(A,b,x0,delta,l)
%用Jacobi迭代法解线线性方程组AX=b,x0为初始向量(以行向量的形式输入),
%delta为给定的误差限,l为给定的最大迭代次数。
[n,m]=size(A);
ifn~=m
fprintf('系数矩阵行列不相等');
return
end
ifn~=length(b)
fprintf('系数矩阵与常数项矩阵的维数不相等');
return
end
fori=1:
n
ifabs(A(i,i))<1e-10
fpritf('失败信息:
因为系数矩阵的主对角线元素近似为0,Jacobi迭代失败。
');
return
end
end
ifnargin<3%如果只输入前两个变量的值,则初值取零向量
x0=zeros(1,n);
end
ifnargin<4
delta=1e-6;
end
ifnargin<5
l=100;
end
fork=1:
l
fori=1:
n
x(i)=(b(i)-A(i,1:
n)*x0'+A(i,i)*x0(i))/A(i,i);
end
ifmax(abs(x-x0)) fprintf('经过%3d次迭代,得到近似解为: ',k); return end x0=x; end fprintf('迭代了%3d次,还是没有达规定的误差要求',l); 2.在MATLAB工作窗口输入: >>A=[12-2;111;221];b=[725]';x0=[000];delta=10^(-5);l=100; >>x=JacobiIteration(A,b,x0,delta,l) 2)运行结果 经过4次迭代,得到近似解为: x= 12-1 G-S迭代法 1)实验编程 1.新建M文件GaussSeidel.m输入: functionx=GaussSeidel(A,b,x0,delta,l) %用高斯-塞德尔迭代法解线线性方程组AX=b,x0为初始向量(以行向量的形式输入), %delta为给定的误差限,l为给定的最大迭代次数。 [n,m]=size(A); ifn~=m fprintf('系数矩阵行列不相等'); return end ifn~=length(b) fprintf('系数矩阵与常数项矩阵的维数不相等'); return end fori=1: n ifA(i,i)==0 fpritf('系数矩阵的主对角元有零元素'); return end end ifnargin<3 x0=zeros(1,n); end ifnargin<4 delta=1e-6; end ifnargin<5 l=100; end x=zeros(1,n); fori=1: n x(i)=(b(i)-A(i,1: i-1)*x(1: i-1)'+A(i,i+1: n)*x0(i+1: n)')/A(i,i); end ifmax(abs(x-x0)) fprintf('经过%3d次迭代,得满足规定精度的近似解为: ',k); x=x0; return end end 2.在MATLAB工作窗口输入: >>A=[311;13-1;31-5];b=[53-1]';x0=[000];delta=10^(-5);l=100; >>x=GaussSeidel(A,b,x0,delta,l) 3)运行结果 经过4次迭代,得到近似解为: x= 1.66670.44441.2889 Euler法 1)实验编程 1.建立名为f.m的M文件输入: functionf=f(x,y) f=1./(x.^2)-y./x; 2.建立名为euler.m的M文件输入: function[x,y]=euler(f,x0,y0,b,n) h=(b-x0)/n; x (1)=x0;y (1)=y0; fork=1: n y(k+1)=y(k)+h*feval(f,x(k),y(k)); x(k+1)=x(k)+h; end 3.在MATLAB工作窗口输入: >>x0=1;y0=1;b=2;n=20; [x,y]=euler('f',x0,y0,b,n) 2)运行结果 x= Columns1through11 1.00001.05001.10001.15001.20001.25001.30001.35001.40001.45001.5000 Columns12through21 1.55001.60001.65001.70001.75001.80001.85001.90001.95002.0000 y= Columns1through11 1.00001.00000.99770.99370.98830.98180.97460.96670.95830.94960.9406 Columns12through21 0.93150.92230.91300.90370.89440.88520.87600.86690.85800.8491 改进Euler法 1)实验编程 1.建立名为f.m的M文件输入: functionf=f(x,y) f=1./(x.^2)-y./x; 2.建立名为EulerModi.m的M文件输入: function[x,y]=EulerModi(f,x0,y0,b,n) h=(b-x0)/n; x (1)=x0;y (1)=y0; fork=1: n x(k+1)=x(k)+h; yy=y(k)+h*feval(f,x(k),y(k)); y(k+1)=y(k)+h/2*(feval(f,x(k),y(k))+feval(f,x(k+1),yy)); end 3.在MATLAB工作窗口输入: >>x0=1;y0=1;b=2;n=20; >>[x,y]=EulerModi('f',x0,y0,b,n) 2)运行结果 x= Columns1through11 1.00001.05001.10001.15001.20001.25001.30001.35001.40001.45001.5000 Columns12through21 1.55001.60001.65001.70001.75001.80001.85001.90001.95002.0000 y= Columns1through11 1.00000.99890.99580.99110.98530.97860.97110.96310.95470.94600.9371 Columns12through21 0.92800.91880.90960.90040.89130.88220.87320.86420.85540.8467 实验结果分析: 1.用Jacobi迭代法可算出方程组的解: (1,2,-1),迭代次数4次。 用G-S迭代法算出方程组的解: (1.6667,0.4444,1.2889),且迭代解比较发散。 2.Euler法和改进的Euler法相比较,改进的Euler法的计算精度更高,相对误差 也比较小。 因此在求解微分方程的数值解时,改进的Euler法优于Euler法。 3.在上述两种方法中当步长h越小则计算精度越高,相对误差较小。 因此,计算能力允许的范围内,选取步长越小可以得到更加精确的结果。 4.在利用这两种方法计算数值解的过程中,前一次迭代的结果会对下一轮求得的数值产生影响。 因此,一旦上一轮迭代所得的结果有偏差,下一轮结果的偏差将大于上一轮的偏差。 因此会导致伴随迭代次数的增加而产生更大的偏差。 成绩评定 签字: 年月日
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值分析 实验报告四 数值 分析 实验 报告