数值分析实验.docx
- 文档编号:4569984
- 上传时间:2023-05-07
- 格式:DOCX
- 页数:53
- 大小:271.30KB
数值分析实验.docx
《数值分析实验.docx》由会员分享,可在线阅读,更多相关《数值分析实验.docx(53页珍藏版)》请在冰点文库上搜索。
数值分析实验
目录
第一章线性方程组的直接解法2
1.1高斯消去法解线性方程组2
1.2列主元高斯消去法解线性方程组3
1.3用LU方法解线性方程组5
1.4利用追赶法求解三对角线性方程组7
第二章解线性方程组的迭代法9
2.1求解线性方程组的Jacobi迭代法9
2.2求解线性方程组的Gauss-Seidel迭代法10
2.3求解线性方程组的SOR迭代法12
第三章解非线性方程(组)的迭代法14
3.1求解非线性方程的一般迭代法14
3.2求解非线性方程的二分法15
3.3求解非线性方程的Newton法17
3.4求解非线性方程组的牛顿法19
第四章多项式插值22
4.1Lagrange插值多项式的构造22
4.2Newton插值多项式的构造23
第五章数值积分27
5.1Simpson公式和梯形公式27
5.2复化Simpson公式和复化梯形公式27
第六章常微分方程初值问题的数值解法29
6.1Euler方法和改进的Euler方法29
6.2四阶龙格-库塔方法27
第一章线性方程组的直接解法
室验内容简介:
如何利用几种常见的直接解法求解线性方程组
实验目的:
加强编程能力和编程技巧,练习从数值分析的角度看问题。
主要仪器及配套软件:
Matlab、C、Mathematica等
1.1高斯消去法解线性方程组
实验项目名称:
用高斯消去法求解下列方程组Ax=b
实验内容:
题目:
已知线性方程组如下:
A=
222
324
139
b=
1.0000
0.5000
2.5000
原理:
将增广矩阵化为下三角阵,从而求解线性方程组。
基本思想:
用Gauss消去法解线性方程组的基本思想是设法消去方程组的系数矩阵A的主对角线下的元素,而将线性方程组化为等价的上三角形方程组,然后再通过回代过程便可获得原方程组的解。
Matlab源代码:
functionx=Gauss_solve(A,b)
%x=Gauss_solve(A,b)returnsthesolutiontotheequationAx=b,
%whereAisann-by-nmatrixandbisacolumnvectorof
%lengthn(oramatrixwithseveralsuchcolumns).
[n,n]=size(A);
fork=1:
n-1
ifA(k,k)==0
disp('Aissingular,stop.')
end
%zerooutentriesofAandbusingpivotA(k,k)
fori=k+1:
n
mik=A(i,k)/A(k,k);
b(i)=b(i)-mik*b(k);
forj=k+1:
n
A(i,j)=A(i,j)-mik*A(k,j);
end
end
end
%backsubstitution
x=zeros(n,1);
x(n)=b(n)/A(n,n);
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A(i,j)*x(j);
end
x(i)=(b(i)-sum)/A(i,i);
end
实验结果:
在命令窗口输入相关数据,运行上面定义的函数Gauss_solve
>>A=[222;324;139]
A=
222
324
139
>>b=[1;1/2;5/2]
b=
1.0000
0.5000
2.5000
>>x=Gauss_solve(A,b)
x=
-0.5000
1.0000
0
1.2列主元高斯消去法解线性方程组
实验项目名称:
用列主元高斯消去法求解下列方程组Ax=b
实验内容:
题目:
A=
0.72900.81000.9000
1.00001.00001.0000
1.33101.21001.1000
b=
0.6867
0.8338
1.0000
原理:
在采用Gauss消去法求解线性方程组时,小主元素可能导致计算失败,故在消去法中应避免采用绝对值很小的主元素。
故在Gauss消去法中的每一步需要采取选主元的技巧。
设计思想:
在高斯消去法的基础上,列主元消去法在每次选主元时,仅依次按列选取绝对值最大的元素作为主元素,且只交换相应的两行(不需要进行列交换),再进行消去计算。
Matlab源代码:
functionx=lsolve(A,b)
%x=lsolve(A,b)returnsthesolutiontotheequationAx=b,
%whereAisann-by-nmatrixandbisacolumnvectorof
%lengthn(oramatrixwithseveralsuchcolumns).
%Gaussianeliminationwithpartialpivoting
[n,n]=size(A);
fork=1:
n-1
%findindexoflargestelementbelowdiagonalincolumnk
ik=k;
fori=k+1:
n
ifabs(A(i,k))>abs(A(ik,k))
ik=i;
end
end
ifA(ik,k)==0
disp('Aissingular,stop.')
end
%swapwithrowk
ifik~=k
forj=k:
n
zjz=A(k,j);
A(k,j)=A(ik,j);
A(ik,j)=zjz;
end
zjz=b(k);
b(k)=b(ik);
b(ik)=zjz;
end
%zerooutentriesofAandbusingpivotA(k,k)
fori=k+1:
n
mik=A(i,k)/A(k,k);
b(i)=b(i)-mik*b(k);
forj=k+1:
n
A(i,j)=A(i,j)-mik*A(k,j);
end
end
end
%backsubstitution
x=zeros(n,1);
x(n)=b(n)/A(n,n);
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A(i,j)*x(j);
end
x(i)=(b(i)-sum)/A(i,i);
end
实验结果:
在命令窗口输入相关数据,运行上面定义的函数lsolve
>>A=[0.7290.8100.900;1.0001.0001.000;1.3311.2101.100]
A=
0.72900.81000.9000
1.00001.00001.0000
1.33101.21001.1000
>>b=[0.6867;0.8338;1.0000]
b=
0.6867
0.8338
1.0000
>>x=lsolve(A,b)
x=
0.2245
0.2814
0.3279
1.3用LU方法解线性方程组
实验项目名称:
用LU法求解下列方程组Ax=b
实验内容:
题目:
A=
222
324
139
b=
1.0000
0.5000
2.5000
原理:
对于线性方程组的系数矩阵A,求出单位下三角矩阵L和上三角矩阵U,使得A=LU。
事实上,L的主对角线左下角元素是求解该线性方程组的Gauss消去过程中使用的乘子,并按消去过程的方式排列所得。
设计思想:
应用LU分解法求解线性方程组Ax=b等价于求解两个三角形线性方程组:
(1)求解下三角线性方程组Ly=b,得y;
(2)求解上三角线性方程组Ux=y,得x。
Matlab源代码:
function[L,U,x]=lu_solve(A,b)
%x=lu_solve(A,b)returnsthesolutiontotheequationAx=b,
%whereAisann-by-nmatrixandbisacolumnvectorof
%lengthn(oramatrixwithseveralsuchcolumns).
%LUdecomposition
%
[n,n]=size(A);
%computethefirstrowofUandthefirstcolumnofL
fori=2:
n
A(i,1)=A(i,1)/A(1,1);
end
fork=2:
n
%computethekthrowofU
forj=k:
n
forp=1:
k-1
A(k,j)=A(k,j)-A(k,p)*A(p,j);
end
end
%computethekthcolumnofL
ifk~n
fori=k+1:
n
forp=1:
k-1
A(i,k)=A(i,k)-A(i,p)*A(p,k);
end
A(i,k)=A(i,k)/A(k,k);
end
end
end
L=zeros(n,n);
fori=1:
n
L(i,i)=1;
end
fori=2:
n
forj=1:
i-1
L(i,j)=A(i,j);
end
end
U=zeros(n,n);
fori=1:
n
forj=i:
n
U(i,j)=A(i,j);
end
end
%solvetheequationsLy=b;
y=b;
fori=2:
n
forj=1:
i-1
y(i)=y(i)-L(i,j)*y(j);
end
end
%solvetheequationsUx=y
x=y;
x(n)=y(n)/U(n,n);
fori=n-1:
-1:
1
forj=i+1:
n
y(i)=y(i)-U(i,j)*x(j);
end
x(i)=y(i)/U(i,i);
end
实验结果:
在命令窗口输入相关数据,运行上面定义的函数lu_solve
>>A=[123;252;315]
A=
123
252
315
>>b=[14;18;20]
b=
14
18
20
>>[L,U,x]=lu_solve(A,b)
L=
100
210
3-51
U=
123
01-4
00-24
x=
1
2
3
1.4利用追赶法求解三对角线性方程组
实验项目名称:
用追赶法求解三对角线性方程组
实验内容:
题目:
用追赶法解方程组
原理:
1.消元过程:
u1=c1/b1,y1=f1/b1
ui=ci/(bi-ui-1*ai)
yi=(fi-yi-1*ai)/(bi-ui-1*ai)
2.回代过程:
xn=yn
xi=yi-ui*xi+1
Matlab源代码:
实验结果:
d=1.0842-0.73680.6000-0.21050.2421
实验注意:
通过本实验要了解追赶法求方程组的解的过程和内涵,这种方法的消元过程和回代过程是关键的地方,通过这次实验要进一步加强对对角占优阵的了解。
第二章解线性方程组的迭代法
实验内容简介:
如何利用几种常见的迭代解法求解线性方程组
实验目的:
练习引入迭代矩阵的形式来解决方程组,巩固方程与矩阵相互关系的概念。
主要仪器及配套软件:
Matlab、C、Mathematica等
2.1求解线性方程组的Jacobi迭代法
实验项目名称:
用Jacobi迭代法求解方程组
实验内容:
题目:
用雅可比迭代求方程组:
AX=b
已知:
A、b如下:
A=
10-120
-111-13
2-110-1
03-18
b=
6
25
-11
15
原理:
先找出下三角阵-L,再找出上对角阵-U,还有主对角阵D。
矩阵形式的迭代公式x=inv(D)*(L+U)*x+inv(D)*b
设计思想:
利用迭代即可在循环中实现的原则来完成此方程组的求解,所用的迭代公式为:
x=inv(D)*(L+U)*x+inv(D)*b
Matlab源代码:
function[iter,x]=Jacobi(A,b,x0,eps,M)
%x0theinitialpoint
%Mthemaximumiteration
%epstheerrortolerance
n=length(x0);
x=x0;
iter=0;
fork=1:
M
fori=1:
n
sum=0;
forj=1:
n
ifj~=i
sum=sum+A(i,j)*x0(j);%x0standsforx_k
end
end
x(i)=(b(i)-sum)/A(i,i);%xstandsforx_{k+1}
end
ifnorm(x-x0) iter=k; return end x0=x; end disp('超过最大迭代次数') 实验结果: 在命令窗口输入相关数据,运行上面定义的函数Jacobi >>A=[10-120;-111-13;2-110-1;03-18] A= 10-120 -111-13 2-110-1 03-18 >>b=[6;25;-11;15] b= 6 25 -11 15 >>x0=[0;0;0;0] x0= 0 0 0 0 >>M=100 M= 100 >>eps=1e-4 eps= 1.0000e-004 >>[iter,x]=Jacobi(A,b,x0,eps,M) iter= 13 x= 1.0000 2.0000 -1.0000 1.0000 2.2求解线性方程组的Gauss-Seidel迭代法 实验项目名称: 用Gauss-Seidel迭代法求解方程组 实验内容: 题目: 用Gauss-Seidel迭代求方程组: AX=b。 方程组如2.1节中的题目。 原理: 先找出下三角阵-L,再找出上对角阵-U,还有主对角阵D。 矩阵形式的迭代公式x=inv(D-L)*U*x+inv(D-L)*b 设计思想: 利用迭代即可在循环中实现的原则来完成此方程组的求解,所用的迭代公式为: x=inv(D-L)U*x+inv(D-L)*b Matlab源代码: function[iter,x]=GS(A,b,x0,eps,M) n=length(x0); iter=0; fork=1: M x=x0;%xstandsforx_{k},x0standsforx_{k+1} fori=1: n ifi==1 sum=0; forj=2: n sum=sum+A(i,j)*x0(j); end t=(b(i)-sum)/A(i,i); x0(i)=t; end%endi==1 if1 sum=0; forj=1: i-1 sum=sum+A(i,j)*x0(j); end forj=i+1: n sum=sum+A(i,j)*x0(j); end t=(b(i)-sum)/A(i,i); x0(i)=t; end%end1 ifi==n sum=0; forj=1: n-1 sum=sum+A(i,j)*x0(j); end t=(b(i)-sum)/A(i,i); x0(i)=t; end%endi==n end%endi ifnorm(x-x0) x=x0; iter=k; return end end%endk disp('超过最大迭代次数') 实验结果: 在命令窗口输入2.1节中一样的相关数据,运行上面定义的函数GS >>[iter,x]=GS(A,b,x0,eps,M) iter= 6 x= 1.0000 2.0000 -1.0000 1.0000 2.3求解线性方程组的SOR迭代法 实验项目名称: 用SOR迭代法求解方程组 实验内容: 题目: 用SOR迭代求方程组: AX=b。 方程组如2.1节中的题目。 原理: 先找出下三角阵-L,再找出上对角阵-U,还有主对角阵D。 矩阵形式的迭代公式x=inv(D-w*L)*[(1-w)*D+w*U]+w*inv(D-w*L)*b,其中w为松弛因子。 设计思想: 利用迭代即可在循环中实现的原则来完成此方程组的求解,所用的迭代公式为: x=inv(D-w*L)*[(1-w)*D+w*U]+w*inv(D-w*L)*b Matlab源代码: function[iter,x]=SOR(w,A,b,x0,eps,M) n=length(x0); iter=0; fork=1: M x=x0; fori=1: n ifi==1 sum=0; forj=2: n sum=sum+A(i,j)*x0(j); end t=(b(i)-sum)/A(i,i); x0(i)=(1-w)*x0(i)+w*t; end%endi==1 if1 sum=0; forj=1: i-1 sum=sum+A(i,j)*x0(j); end forj=i+1: n sum=sum+A(i,j)*x0(j); end t=(b(i)-sum)/A(i,i); x0(i)=(1-w)*x0(i)+w*t; end%end1 ifi==n sum=0; forj=1: n-1 sum=sum+A(i,j)*x0(j); end t=(b(i)-sum)/A(i,i); x0(i)=(1-w)*x0(i)+w*t; end%endi==n end%endi ifnorm(x0-x) x=x0; iter=k; return end end%endk disp('超过最大迭代次数') 实验结果: 在命令窗口输入2.1节中一样的相关数据,运行上面定义的函数SOR >>w=1.2 w= 1.2000 >>[iter,x]=SOR(w,A,b,x0,eps,M) iter= 7 x= 1.0000 2.0000 -1.0000 1.0000 >>w=1.5 w= 1.5000 >>[iter,x]=SOR(w,A,b,x0,eps,M) iter= 15 x= 1.0000 2.0000 -1.0000 1.0000 第三章解非线性方程(组)的迭代法 实验内容简介: 如何利用几种常见的迭代解法求解非线性方程(组)。 实验目的: 通过本章的练习,加深理解迭代法的真谛,掌握其中的规律。 主要仪器及配套软件: Matlab、C、Mathematica等 3.1求解非线性方程的一般迭代法 实验项目名称: 通过构造简单迭代格式求解非线性方程 实验内容: 题目: 给出计算x= 的迭代公式,讨论迭代过程的收敛性,并证明x=2。 原理: 根据上式的规律可以看出,若将每个根号内的表达式看作一个整体,则这个整体就会呈现出迭代的规律。 设计思想: 经以上分析,令x=(2+x)^0.5作为迭代公式。 Matlab源代码: disp('给出计算x=(2+(2+(2+((2+...)^0.5)^0.5)^0.5)^0.5)^0.5') disp('的迭代公式,讨论迭代过程的收敛性,并证明x=2') fprintf('\n迭代公式为x=(2+x)^0.5\n') x1=2^0.5; x2=2; fazi=input('请输入精度(当精度为0时,说明为精确值! )>>') whileabs(x2-x1)>fazi x1=(2+x1)^0.5; x2=(2+x1)^0.5; end ifx1==2 fprintf('\n此时为准确值: x=%f\n',x1) else fprintf('\n近似值: x=%f',x1) end fprintf('\n/*******************************************/') 实验结果: 3.2求解非线性方程的二分法 实验项目名称: 利用二分法求解非线性方程 实验内容: 题目: 在区间[1,3]用二分法求 。 原理: 将给定的初始有根区间逐步分半,直到找到满足一定精度要求的解。 设计思想: 对当前的有根区间,取其中点作为真解的近似值,若满足精度要求则停止计算,否则,从得到的两个小区间中选择一个作为新的有根区间,重复上述过程。 程序中重要的是循环语句的使用和有根区间端点的判断。 Matlab源代码: functionf=fun(x) f=2^(-x)+exp(x)+2*cos(x)-6; function[iter,x]=bisection(a,b,eps) %abarethegiveninitialpointssuchthatf(a)*f(b)<0 %epsistheerrortolerance funa=fun(a); funb=fun(b); iffuna*funb>=0 disp('Error,pleasegivethecorrectinitialpointsaandb! '); end k=0;%kstandsfortheiteration while1 c=(a+b)/2; func=fun(c); k=k+1; if(abs(c-a)<=eps)|(abs(func)<=eps) x=c; iter=k; return; end iffuna*func<0 b=c; else
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验