《数值分析》课程实验报告

　 姓 名：

　学 号：

　学 院：

　机 电 学 院

　日 期：

　 2015 年 X 月X 日

　目

　录 实验一 函数插值方法 1 实验二 函数逼近与曲线拟合 5 实验三 数值积分与数值微分 7 实验四 线方程组的直接解法 9 实验五 解线性方程组的迭代法 15 实验六 非线性方程求根 19 实验七 矩阵特征值问题计算 21 实验八 常微分方程初值问题数值解法 24

　 实验一 函数插值方法 一、问题提出

　对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。

　数据如下：

　 （1）

　0.4

　0.55

　0.65

　0.80

　0.95

　1.05

　 0.41075

　0.57815 0.69675 0.90

　1.00

　1.25382

　 求五次Lagrange多项式，和分段三次插值多项式，计算, 的值。（提示：结果为,

　 ）

　（2）

　1

　2

　3

　4

　5

　6

　7

　 0.368

　0.135

　0.050

　0.018

　0.007

　0.002

　0.001

　试构造Lagrange多项式，计算的,值。（提示：结果为,

　 ） 二、要求

　1、 利用Lagrange插值公式

　 编写出插值多项式程序；

　2、 给出插值多项式或分段三次插值多项式的表达式；

　3、 根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何；

　4、 对此插值问题用Newton插值多项式其结果如何。Newton插值多项式如下：

　其中：

　 三、目的和意义

　1、 学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；

　2、 明确插值多项式和分段插值多项式各自的优缺点；

　3、 熟悉插值方法的程序编制；

　4、 如果绘出插值函数的曲线，观察其光滑性。

　四、实验步骤 （1）

　0.4

　0.55

　0.65

　0.80

　0.95

　1.05

　 0.41075

　0.57815 0.69675 0.90

　1.00

　1.25382

　 求五次Lagrange多项式，和分段三次插值多项式，计算, 的值。（提示：结果为,

　 ）

　第一步：先在matlab中定义lagran的M文件为拉格朗日函数代码为：

　 function[c,l]=lagran(x,y) w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if(k~=j) v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; end

　第二步：然后在matlab命令窗口输入：

　>>>> x=[0.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382]; >> lagran(x,y)

　回车得到：

　 ans =121.6264

　-422.7503

　572.5667 -377.2549

　121.9718

　-15.0845 由此得出所求拉格朗日多项式为

　p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845

　第三步：在编辑窗口输入如下命令：

　>> x=[0.4 0.55 0.65 0.80,0.95 1.05]; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.0845; >> plot(x,y)

　命令执行后得到如下图所示图形，然后 >> x=0.596; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084 y =0.6262

　 得到f（0.596）=0.6262

　 同理得到f（0.99）=1.0547

　 （2）

　1

　 2

　3

　4

　5

　6

　7

　 0.368

　 0.135

　0.050

　0.018

　0.007

　0.002

　0.001

　试构造Lagrange多项式，和分段三次插值多项式，计算的,值。（提示：结果为,

　 ）

　 实验步骤：

　第一步定义 function[c,l]=lagran(x,y) w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if(k~=j) v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; end 定义完拉格朗日M文件

　 第二步：然后在matlab命令窗口输入：

　>>>> x=[1 2 3 4 5 6 7]; y=[0.368 0.135 0.050 0.018 0.007 0.002 0.001]; >> lagran(x,y) 回车得到：

　ans =0.0001

　 -0.0016

　0.0186

　 -0.1175

　0.4419

　 -0.9683

　0.9950 由此得出所求拉格朗日多项式为 p（x）=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950 第三步：在编辑窗口输入如下命令：

　 >> x=[1 2 3 4 5 6 7]; >> y=0.0001*x.^6-0.0016*x.^5+0.0186*x.^4-0.1175*x.^3+0.4419*x.^2-0.9683*x+0.9950; >> plot(x,y) 命令执行后得到如下图所示图形，然后 >> x=1.8; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084 y =0.1650

　 得到f（0.596）=0.6262

　同理得到f（6.15）=2.3644

　 五、实验结论

　插值是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点，它是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。

　 实验二 函数逼近与曲线拟合

　一、问题提出

　从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。

　 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。

　t(分) 0

　 5

　10

　15

　20

　25

　30

　35

　40

　45

　50

　55

　 0

　1.27

　2.16

　2.86

　3.44

　3.87

　4.15

　4.37

　4.51

　4.58

　4.02

　4.64

　 二、要求

　1、用最小二乘法进行曲线拟合；

　2、近似解析表达式为； 3、打印出拟合函数，并打印出与的误差，；

　4、另外选取一个近似表达式，尝试拟合效果的比较；

　5、* 绘制出曲线拟合图。

　三、目的和意义

　1、掌握曲线拟合的最小二乘法；

　2、最小二乘法亦可用于解超定线代数方程组；

　3、探索拟合函数的选择与拟合精度间的关系

　四、实验步骤：

　 第一步先写出线性最小二乘法的M文件 function c=lspoly(x,y,m) n=length(x); b=zeros(1:m+1); f=zeros(n,m+1); for k=1:m+1 f(:,k)=x.^(k-1); end a=f＇*f; b=f＇*y＇; c=a＼b; c=flipud(c);

　第二步在命令窗口输入：

　>>lspoly([0,5,10,15,20,25,30,35,40,45,50,55],[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64],2)

　回车得到：

　ans =

　 -0.0024

　 0.2037 0.2305 即所求的拟合曲线为y=-0.0024x2+0.2037x+0.2305 在编辑窗口输入如下命令：

　>> x=[0,5,10,15,20,25,30,35,40,45,50,55]; >> y=-0.0024*x.^2+0.2037*x+0.2305; >> plot(x,y) 命令执行得到如下图

　 五、实验结论  分析复杂实验数据时,常采用分段曲线拟合方法。利用此方法在段内可以实现最佳逼近,但在段边界上却可能不满足连续性和可导性。分段函数的光滑算法,给出了相应的误差分析.给出了该方法在分段曲线拟合中的应用方法以及凸轮实验数据自动分段拟合。

　 实验三 数值积分与数值微分

　一、问题提出

　选用复合梯形公式，复合Simpson公式，Romberg算法，计算

　（1）

　（2）

　（3）

　（4）

　 二、要求

　1、 编制数值积分算法的程序；

　2、 分别用两种算法计算同一个积分，并比较其结果；

　3、 分别取不同步长，试比较计算结果（如n = 10, 20等）；

　4、 给定精度要求ε，试用变步长算法，确定最佳步长。

　三、目的和意义

　1、 深刻认识数值积分法的意义；

　2、 明确数值积分精度与步长的关系；

　3、 根据定积分的计算方法，可以考虑二重积分的计算问题。

　 四、 实验步骤

　第一步：编写各种积分的程序

　复合梯形程序如下：

　 function I=TX(x,y)

　n=length(x);m=length(y);

　if n~=m

　 error(＇The lengths of X and Y must be equal＇);

　return;

　end

　h=(x(n)-x(1))/(n-1);

　a=[1 2*ones(1,n-2) 1];

　I=h/2*sum(a.*y);

　 复合Simpson程序如下：

　 function s = simpr1(f,a,b,n)

　h=(b-a)/(2*n);

　s1=0;

　s2=0;

　for k=1:10

　x=a+h*(2*k-1);

　s1=s1+feval(f,x);

　end

　for k=1:(10-1)

　x=a+h*2*k;

　s2=s2+feval(f,x);

　end

　s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;

　end

　 Romberg程序如下：

　function I = Romber_yang(fun,a,b,ep) if nargin<4

　 ep=1e-5;end; m=1;

　h=b-a; I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I; while 1 N=2^(m-1);h=h/2;I=I/2; for i=1:N I=I+h*feval(fun,a+(2*i-1)*h); end T(m+1,1)=I;M=2*N;k=1; while M>1; T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1); M=M/2;k=k+1; end if abs(T(k,k)-T(k-1,k-1))<ep break; end m=m+1; end I=T(k,k);

　第二步：用复合梯形公式，复合Simpson公式，Romberg公式对各个积分进行计算 1、 对于积分Ι=0144-sin2xdx，梯形积分T=0.49871011，辛普森积分S=0.49871111，Romberg积分R=0.49871111。

　2、 对于积分Ι=01sin⁡XXdx，f(0)=1,梯形积分T=0.94607307，辛普森积分S=0.94607308，Romberg积分R=0.94607307。

　3、 对于积分Ι=01eX4+X2dx，梯形积分T=0.39081248，辛普森积分S=0.39081185，Romberg积分R=0.39081885。

　4、 对于积分Ι=01ln1+X1+X2dx，梯形积分T=0.27218912，辛普森积分S=0.27219844，Romberg积分R=0.27219827。

　 五、实验结论 ， 通过本实验学会复合梯形公式，复合Simpson公式，Romberg公式的编程与应用，掌握MATLAB提供的计算积分的各种函数的使用方法。

　实验四 线方程组的直接解法 一、问题提出

　给出下列几个不同类型的线性方程组，请用适当算法计算其解。

　 1、 设线性方程组

　 

　2、 设对称正定阵系数阵线方程组

　

　 3、 三对角形线性方程组

　 二、要求

　1、 对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）；

　2、 应用结构程序设计编出通用程序；

　3、 比较计算结果，分析数值解误差的原因；

　4、 尽可能利用相应模块输出系数矩阵的三角分解式。

　三、目的和意义

　1、通过该课题的实验，体会模块化结构程序设计方法的优点；

　2、运用所学的计算方法，解决各类线性方程组的直接算法；

　3、提高分析和解决问题的能力，做到学以致用；

　4、 通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。

　四、实验步骤：

　 列主元高斯消去法的matlab的M文件程序 function [x,det,index]=Gauss(A,b) % 求线形方程组的列主元Gauss消去法，其中， % A为方程组的系数矩阵； % b为方程组的右端项； % x为方程组的解； % det为系数矩阵A的行列式的值； % index为指标变量，index=0表示计算失败，index=1表示计算成功。

　[n,m]=size(A); nb=length(b); % 当方程组行与列的维数不相等时，停止计算，并输出出错信息。

　if n~=m

　 error(＇The rows and columns of matrix A must be equal!＇);

　 return; end % 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息 if m~=nb

　 error(＇The columns of A must be equal the length of b!＇);

　 return; end % 开始计算，先赋初值 index=1;det=1;x=zeros(n,1); for k=1:n-1

　 % 选主元

　 a_max=0;

　 for i=k:n

　 if abs(A(i,k))>a_max

　 a_max=abs(A(i,k));r=i;

　 end

　 end

　 if a_max<1e-10

　 index=0;return;

　 end

　 % 交换两行

　 if r>k

　 for j=k:n

　 z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;

　 end

　 z=b(k);b(k)=b(r);b(r)=z;det=-det;

　 end

　 % 消元过程

　 for i=k+1:n

　 m=A(i,k)/A(k,k);

　 for j=k+1:n

　 A(i,j)=A(i,j)-m*A(k,j);

　 end

　 b(i)=b(i)-m*b(k);

　 end

　 det=det*A(k,k); end det=det*A(n,n); % 回代过程 if abs(A(n,n))<1e-10

　 index=0;return; end for k=n:-1:1

　 for j=k+1:n

　 b(k)=b(k)-A(k,j)*x(j);

　 end

　 x(k)=b(k)/A(k,k); end 然后在命令窗口输入 >> A=[4 2 -3 -1 2 1 0 0 0 0;8 6 -5 -3 6 5 0 1 0 0;4 2 -2 -1 3 2 -1 0 3 1;0 -2 1 5 -1 3 -1 1 9 4;-4 2 6 -1 6 7 -3 3 2 3;8 6 -8 5 7 17 2 6 -3 5;0 2 -1 3 -4 2 5 3 0 1;16 10 -11 -9 17 34 2 -1 2 2;4 6 2 -7 13 9 2 0 12 4;0 0 -1 8 -3 -24 -8 6 3 -1]; >> b=[5 12 3 2 3 46 13 38 19 -21]; >> gauss(A,b) ans =

　1.0000

　-1.0000

　 0.0000

　 1.0000

　 2.0000

　 0.0000

　 3.0000

　 1.0000

　-1.0000 2.0000

　 高斯-约当消去法maltab的M文件程序 function [x,flag]=Gau_Jor(A,b) % 求线形方程组的列主元Gauss-约当法消去法，其中， % A为方程组的系数矩阵； % b为方程组的右端项； % x为方程组的解； [n,m]=size(A); nb=length(b); % 当方程组行与列的维数不相等时，停止计算，并输出出错信息。

　if n~=m

　 error(＇The rows and columns of matrix A must be equal!＇);

　 return; end % 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息 if m~=nb

　 error(＇The columns of A must be equal the length of b!＇);

　 return; end % 开始计算，先赋初值 flag=＇ok＇;x=zeros(n,1); for k=1:n

　 % 选主元

　 max1=0;

　 for i=k:n

　 if abs(A(i,k))>max1

　 max1=abs(A(i,k));r=i;

　 end

　 end

　 if max1<1e-10

　 falg=＇failure＇;return;

　 end

　 % 交换两行

　 if r>k

　 for j=k:n

　 z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;

　 end

　 z=b(k);b(k)=b(r);b(r)=z;

　 end

　 % 消元过程

　 b(k)=b(k)/A(k,k);

　 for j=k+1:n

　 A(k,j)=A(k,j)/A(k,k);

　 end

　 for i=1:n

　 if i~=k

　 for j=k+1:n

　 A(i,j)=A(i,j)-A(i,k)*A(k,j);

　 end

　 b(i)=b(i)-A(i,k)*b(k);

　 end

　 end end % 输出x for i=1:n

　 x(i)=b(i); end 然后保存后在命令窗口输入：

　>> A=[4 2 -4 0 2 4 0 0;2 2 -1 -2 1 3 2 0;-4 -1 14 1 -8 -3 5 6;0 -2 1 6 -1 -4 -3 3;2 1 -8 -1 22 4 -10 -3;4 3 -3 -4 4 11 1 -4;0 2 5 -3 -10 1 14 2;0 0 6 3 -3 -4 2 19]; >> b=[0 -6 20 23 9 -22 -15 45]; >> Gau_Jor(A,b) ans =

　121.1481

　-140.1127

　29.7515

　 -60.1528

　10.9120

　 -26.7963

　 5.4259

　-2.0185

　 五、实验结论 用LU法，调用matlab中的函数lu中，L往往不是一个下三角，但可以直接计算不用它的结果来计算，不用进行行变换。如果进行行变b也要变，这样会很麻烦。

　 实验五 解线性方程组的迭代法 一、问题提出

　对实验四所列目的和意义的线性方程组，试分别选用Jacobi 迭代法，Gauss-Seidel迭代法和SOR方法计算其解。

　二、要求

　1、体会迭代法求解线性方程组，并能与消去法做以比较；

　2、分别对不同精度要求，如由迭代次数体会该迭代法的收敛快慢； 3、对方程组2，3使用SOR方法时，选取松弛因子ω=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者；

　4、给出各种算法的设计程序和计算结果。

　三、目的和意义

　1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；

　2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序；

　3、体会上机计算时，终止步骤或k >（给予的迭代次数），对迭代法敛散性的意义；

　4、 体会初始解，松弛因子的选取，对计算结果的影响。

　 四、实验步骤 第一步编写实验所需的Jacobi迭代法，Gauss-Seidel迭代法，SOR迭代法的程序。

　Jacobi迭代法：

　function [x,k,index]=J(A,b,ep,itmax) if nargin<4

　 itmax=100; end if nargin<3 ep=1e-5; end n=length(A); k=0; x=zeros(n,1);y=zeros(n,1);index=1; while 1

　 for i=1:n

　 y(i)=b(i);

　 for j=1:n

　 if j~=i

　 y(i)=y(i)-A(i,j)*x(j);

　 end

　 end

　 if abs(A(i,i))<1e-10|k==itmax

　 index=0;return; end y(i)=y(i)/A(i,i); end if norm(y-x,inf)<ep

　 break; end x=y;k=k+1; end

　Gauss-Seidel迭代法：

　function [x,k,index]=G(A,b,ep,itmax) if nargin<4

　 itmax=100; end if nargin<3 ep=1e-5; end n=length(A); k=0; x=zeros(n,1);y=zeros(n,1);index=1; while 1

　 y=x;

　 for i=1:n

　 z=b(i);

　 for j=1:n

　 if j~=i

　 z=z-A(i,j)*x(j);

　 end

　 end

　 if abs(A(i,i))<1e-10|k==itmax

　 index=0;return;

　 end

　 z=z/A(i,i);x(i)=z;

　 end

　 if norm(y-x,inf)<ep

　 break;

　 end

　 k=k+1; end

　SOR迭代法：

　function [x,k,index]=SOR(A,b,ep,w,itmax) if nargin<5

　 itmax=100; end if nargin<4 w=1;end if nargin<3 ep=1e-5; end n=length(A); k=0; x=zeros(n,1);y=zeros(n,1);index=1; while 1

　 y=x;

　 for i=1:n

　 z=b(i);

　 for j=1:n

　 if j~=i

　 z=z-A(i,j)*x(j);

　 end

　 end

　 if abs(A(i,i))<1e-10|k==itmax

　 index=0;return;

　 end

　 z=z/A(i,i);x(i)=(1-w)*x(i)+w*z;

　 end

　 if norm(y-x,inf)<ep

　 break;

　 end

　 k=k+1; end

　第二步：把实验给出的方程代入程序计算其结果。

　 1、 设线性方程组

　 

　2、 设对称正定阵系数阵线方程组

　3、三对角形线性方程组

　 五、实验结论 迭代法是解线性方程组的一个重要的实用方法，特别适用于求解在实际中大量出现的，系数矩阵为稀疏阵的大型线性方程组。通过此次实验学会了Jacobi迭代法，Gauss-Seidel迭代法，SOR迭代法的程序编写，并掌握了它们各自的优缺点及其适用条件。

　 实验六 非线性方程求根 一、问题提出

　设方程有三个实根 现采用下面六种不同计算格式，求 f(x)=0的根或 1、

　2、

　3、

　4、

　5、

　6、

　 二、要求

　1、编制一个程序进行运算，最后打印出每种迭代格式的敛散情况；

　2、用事后误差估计来控制迭代次数，并且打印出迭代的次数；

　3、初始值的选取对迭代收敛有何影响；

　4、分析迭代收敛和发散的原因。

　三、目的和意义

　1、通过实验进一步了解方程求根的算法；

　2、认识选择计算格式的重要性；

　3、掌握迭代算法和精度控制；

　4、明确迭代收敛性与初值选取的关系。

　四、实验步骤 第一步：编写实验所需的程序。

　function [x_star,index,it]=DD(fun,x,ep,itmax) if nargin<4

　 itmax=100; end if nargin<3 ep=1e-5; end index=0;k=1; while k<itmax

　 x1=x;x=feval(fun,x);

　 if abs(x-x1)<ep

　 index=1;break;

　 end

　 k=k+1; end x_star=x;it=k;

　第二步：分别用上述六种形式的表达式计算方程的根，结果如下。

　1、 ，x1=0，x2=0。

　2、 ，x1=无穷大，x2=-0.3473。

　3、 ，x1=1.8794，x2=1.8794。

　4、 ，x1=-0.3473，x2=-0.3473.。

　5、 ，x1=1.8794，x2=1.8794。

　6、 ，x1=1.8794，x2=-0.3473。

　 五、 实验结论 对于非线性方程，求它的解析解有时候是很困难的，但采用数值方法可以很容易地求它的近似解。此次实验就是采用迭代法求非线性方程的根。对于一个非线性方程，选用不同的迭代形式，因为其收敛程度不一样，造成其效率与精确度有很大的差别。

　实验七 矩阵特征值问题计算 一、问题提出

　利用冪法或反冪法，求方阵的按模最大或按模最小特征值及其对应的特征向量。

　 设矩阵A的特征分布为：

　 且

　试求下列矩阵之一

　（1） 求，及

　 取

　 结果 （2） 求及

　 取

　 结果：

　（3） 求及 取 结果

　(4)

　 取

　 这是一个收敛很慢的例子,迭代次才达到 结果

　（5）

　 有一个近似特征值,试用幂法求对应的特征向量，并改进特征值（原点平移法）。

　取

　 结果

　二、要求

　1、掌握冪法或反冪法求矩阵部分特征值的算法与程序设计；

　2、会用原点平移法改进算法，加速收敛；对矩阵B=A-PI取不同的P值，试求其效果；

　3、试取不同的初始向量，观察对结果的影响；

　4、对矩阵特征值的其它分布，如且如何计算。

　三、目的和意义

　1、求矩阵的部分特征值问题具有重要实际意义，如求矩阵谱半径，稳定性问题往往归于求矩阵按模最小特征值；

　2、进一步掌握冪法、反冪法及原点平移加速法的程序设计技巧；

　3、问题中的题（5），反应了利用原点平移的反冪法可求矩阵的任何特征值及其特征向量。

　四、实验步骤

　第一步：写出实验所需的幂法求最大特征值及反幂法求最小特征值的程序。

　幂法程序：

　function [m,u,index]=TZ(A,ep,itmax) if nargin<3

　 itmax=100; end if nargin<2 ep=1e-5; end n=length(A); u=ones(n,1); index=0;k=0;m1=0; while k<=itmax

　 v=A*u;[vmax,i]=max(abs(v));

　 m=v(i);u=v/m;

　 if abs(m-m1)<ep

　 index=1;break;

　 end

　 m1=m;k=k+1; end

　反幂法程序：

　function [m,u,index]=FTZ(A,ep,itmax) if nargin<3

　 itmax=100; end if nargin<2 ep=1e-5; end n=length(A); u=ones(n,1); index=0;k=0;m1=0; invA=inv(A); while k<=itmax

　 v=invA*u;[vmax,i]=max(abs(v));

　 m=v(i);u=v/m;

　 if abs(m-m1)<ep

　 index=1;break;

　 end

　 m1=m;k=k+1; end

　 第二步：将所给矩阵代入程序得出结果。

　，λ1=-6.4210，X1=（-0.0462 -0.3749 1.0000）T；λ3=3.4723，x3=（1.0000 0.5229 0.2422）T。

　，λ1=21.3053，X1=（0.8724 0.5401 0.9973 0.5644 0.4972 1.0000）T；λ6=1.6214。

　 五、实验结论 求n阶方阵A的特征值和特征向量，也是实际中常常碰到的问题。通过此次实验掌握了用幂法和反幂法求一个方阵的最大特征值和特征向量，绝对值最小的特征值和特征向量。

　 实验八 常微分方程初值问题数值解法 一、问题提出

　科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法，Rung-Kutta方法求其数值解，诸如以下问题：

　 （1）

　分别取h=0.1,0.2,0.4时数值解。

　初值问题的精确解。

　 （2）

　用r=3的Adams显式和预 - 校式求解

　 取步长h=0.1，用四阶标准R-K方法求值。

　 （3）

　用改进Euler法或四阶标准R-K方法求解

　取步长0.01,计算数值解，参考结果 。

　 （4）利用四阶标准R- K方法求二阶方程初值问题的数值解

　（I）

　(II)

　(III)

　 (IV) 

　 二、要求

　1、 根据初值问题数值算法，分别选择二个初值问题编程计算；

　2、 试分别取不同步长，考察某节点处数值解的误差变化情况；

　3、 试用不同算法求解某初值问题，结果有何异常；

　4、 分析各个算法的优缺点。

　三、目的和意义

　1、 熟悉各种初值问题的算法，编出算法程序；

　2、 明确各种算法的精度与所选步长有密切关系；

　3、 通过计算更加了解各种算法的优越性。

　 四、实验步骤

　function [x,y]=euler(fun,x0,xfinal,y0,n); if nargin<5,n=50; end h=(xfinal-x0)/n; x(1)=x0;y(1)=y0; for i=1:n;

　 x(i+1)=x(i)+h;

　 y(i+1)=y(i)+h*feval(fun,x(i),y(i)); end 实验程序及分析

　（Ⅰ）

　（1）、算法程序

　function E =Euler_1(fun,x0,y0,xN,N)

　 % Euler向前公式，其中

　% fun为一阶微分方程的函数

　 % x0,y0为初始条件

　% xN为取值范围的一个端点

　% h为区间步长

　% N为区间个数

　% x为Xn构成的向量

　 % y为yn构成的向量

　x=zeros(1,N+1);y=zeros(1,N+1);

　 x(1)=x0;y(1)=y0;

　 h=(xN-x0)/N;

　 for n=1:N

　x(n+1)=x(n)+h;

　y(n+1)=y(n)+h*feval(fun,x(n),y(n));

　 end

　T=[x＇,y＇]

　function z=f(x,y)

　 z=4*x/y-x*y;

　(2)、运行程序

　 >> Euler_1(＇f＇,0,3,2,20)

　 结果 ：

　 >>

　Euler_1(＇f＇,0,3,2,20)

　T =

　0

　3.0000

　0.1000

　2.9836

　0.2000

　2.9517

　0.3000

　2.9058

　0.4000

　2.8481

　0.5000

　2.7810

　0.6000

　2.7073

　0.7000

　2.6297

　0.8000

　2.5511

　0.9000

　2.4739

　1.0000

　2.4004

　1.1000

　2.3325

　1.2000

　2.2714

　1.3000

　2.2177

　1.4000

　2.1717

　1.5000

　2.1332

　1.6000

　2.1017

　1.7000

　2.0765

　1.8000

　2.0567

　1.9000

　2.0414

　2.0000

　2.0299

　 五、实验结论 很多科学技术和工程问题常用微分方程的形式建立数学模型，因此微分方程的求解是很有意义的，但对于绝大多数的微分方程问题很难或者根本不可能得到它的解析解，因此，研究微分方程的数值方法是非常有意义的。