数值分析大作业QR分解

题目:设计求取n n ?实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵

20010-1-24A=0-2131

43

1?? ? ? ? ???

为例进行具体的求解计算。

一、 算法分析：

一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。

二、 算法设计：

1、 原始实矩阵A R

n n

?∈拟上三角化

为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’

记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)r

ij a i n j r r n ==+。 对于r=1,2,…,n -2执行

(1) 若()

(2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。

(2) 计算 令

()2

()()

1,1,1,sgn(0,sgn()=1)

r r r r r r

r r r r r c a

a a ρ+++=-=，（若则取

(3) 令-0=r n r u u ?? ?

??

，()()()-1,2,1(,,...,)r r r T

n r r r r r r nr u a c a a ρ++=- (4) 计算

(r)(r)(r)T

n-r r+1,r r r+2,r nr r T

n-r T n-r n-r n-r n-r r+1r 1u =

(a -c ,a ,...,a )ρ

I H =I -2uu =H H =I -2u u A =HA H

??

?

??

(5) 继续

算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。 2、 拟上三角矩阵QR 分解的求原矩阵的全部特征值

记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵，

对于k =1,2,…,n -1执行： (1) 分解k-1k-1k-1A =R Q

选取旋转矩阵1P =R(2,1,θ)，使(1)1k-1A =PA ，从而使第一列次对角元(1)

2,1=0a ；再选取旋转矩阵2P =R(3,2,θ)，使(2)(1)1

A =PA ，从而使第一列次对角元(2)3,2=0a ……如此进行下去，最多经过n-1步，(n-1)

A

必然变为上三角矩阵k-1R ，即

k-1n-121k-1

-1=

-1-1

-1

k-1n-12112

n-1

R =P P P A =P P P P P P

Q （）

k-1Q 必为正交矩阵，且满足k-1k-1k-1A =R Q

(2) 计算k k-1k-1A =R Q

(3) 上述两过程反复进行直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。 3、 带位移的QR 方法求实矩阵A R

n n

?∈全部特征值

一般情况下，QR 分解的收敛速度比较慢，因此可通过仿乘幂法将原矩阵先进行一定长度的位移，可显著加速QR 算法所得矩阵序列k A 的收敛。0A =A （A 是原始矩阵的拟上三角矩阵）

(1) 分解k-1-1k-1k-1A -u I=R k Q ，其中-1u k 即位移量，一般取A 的某一特征值的近似值；

实际计算通常取为k-1A 的右下角元素(k-1)n,n a ，或取为右下角

22?矩阵特征之中接近(k-1)n,n a 者。

(2) 计算k k-1k-1-1A =+R u I k Q 。

(3) 重复过程（1）（2）直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。 4、 反幂法求实矩阵A R

n n

?∈的特征向量

记通过QR 分解得到A 的某一特征值的近似值为p ，反幂法步骤如下： (1) 三角分解A -pI =LR 。 (2) 对k=1，2，…，做

①当k=1时，令T

u =(1 11) 当k 1≠时，解-1Lu=z k

②解Ry =u k

③求-1y k 绝对值最大的元素k m ④令k k k z =y /m

⑤当k k-1m -m 或k k-1z -z 小于规定的误差界时，停止计算，k z 即为所求的特征向量，k p +1/m 即为对应特征值的更为精确的取值。

三、 程序设计：

1、 对生成的矩阵A 进行拟上三角化

利用hohd 函数对A 进行householder 变换，得到A 的拟上三角矩阵。

2、 对拟上三角化后的矩阵进行QR 分解和带位移的QR 分解，求矩阵的全部特征值

在绝对误差界为

-41

102

?的条件下，利用qrfenjie 函数对拟上三角矩阵进行QR 分解，利用dp_qrfenjie 函数进行带位移的QR 分解，比较两者的收敛速度。 3、 反幂法求更精确的特征值和特征向量

利用vect 函数和（2）得到的特征值的近似值计算更为精确的特征值和对应的特征向量，是绝对误差界为

-61

102

?。 4、 输出相关结果。

四、 源程序：

1、 hohd 函数

function[A]=hohd(a) n=length(a); for i=1:n-2

b=a(i+1:n,i); if b(1)>=0

c=-sqrt(sum(b.^2)); else

c=sqrt(sum(b.^2)); end

rho=sqrt(2*c*(c-b(1))); u1=(b-c*eye(n-i,1))/rho; H1=eye(n-i)-2*u1*u1'; H=eye(n-i);

H(i+1:n,i+1:n)=H1; a=H*a*(H'); end

A=a

2、qrfenjie函数

function A=qrfenjie(a)

n=length(a);

A=a;

for i=1:100

theta(n-1)=0;

c(n-1)=0;

s(n-1)=0;

Q=1;

for j=1:n-1

theta(j)=atan(A(j+1,j)/A(j,j));

c(j)=cos(theta(j));

s(j)=sin(theta(j));

P=eye(n);

P(j,j)=c(j);

P(j+1,j)=-s(j);

P(j,j+1)=s(j);

P(j+1,j+1)=c(j);

invP=eye(n);

invP(j,j)=c(j);

invP(j+1,j)=s(j);

invP(j,j+1)=-s(j);

invP(j+1,j+1)=c(j);

Q=Q*invP;

R=P*A;

A=R;

end

A=R*Q;

tor=max(abs(diag(A)-diag(a)));

if tor>5*10^-5

a=A;

else

break;

end

end

i,Ak=A,lanmda=diag(A)'

3、dp_qrfenjie函数

function A=dp_qrfenjie(a)

n=length(a);

A=a;

A=a-a(n,n)*eye(n);

theta(n-1)=0;

c(n-1)=0;

s(n-1)=0;

Q=1;

for j=1:n-1

theta(j)=atan(A(j+1,j)/A(j,j));

c(j)=cos(theta(j));

s(j)=sin(theta(j));

P=eye(n);

P(j,j)=c(j);

P(j+1,j)=-s(j);

P(j,j+1)=s(j);

P(j+1,j+1)=c(j);

invP=eye(n);

invP(j,j)=c(j);

invP(j+1,j)=s(j);

invP(j,j+1)=-s(j);

invP(j+1,j+1)=c(j);

Q=Q*invP;

R=P*A;

A=R;

end

A=R*Q+a(n,n)*eye(n);

tor=max(abs(diag(A)-diag(a)));

if tor>5*10^-5

a=A;

else

break;

end

end

i,Ak=A,lanmda=diag(A)'

4、vect函数

function z=vect(a0,p)

n=length(a0);

Z=zeros(n,n+2);

for x=1:n

a=a0-p(x)*eye(n);

r=zeros(n,n);

r(1,1:n)=a(1,1:n);

l(2:n,1)=a(2:n,1)/a(1,1);

for i=2:n

for j=i:n

M=0;

for k=1:i-1

M=l(i,k)*r(k,j)+M;

end

r(i,j)=a(i,j)-M;

end

if i

for j=i+1:n

N=0;

for k=1:i-1

N=l(j,k)*r(k,i)+N;

end

l(j,i)=(a(j,i)-N)/r(i,i);

end

else

break;

end

end

R=r;

L=l;

u=ones(n,1);

mo=0;

invr=inv(R);

invl=inv(L);

for i=1:100

y=invr*u;

p1=max(y);

p2=max(-y);

if p1>p2

m=p1;

else

m=-p2;

end

z=y/m;

tor=abs(m-mo);

if tor<5*10^-7

break;

else

mo=m;

end

u=invl*z;

end

Z(x,3:n+2)=z';

Z(x,1)=i;

Z(x,2)=p(x)+1/m;

end

Z

五、运行结果：

1、执行命令：

>> clear,clc,a=[2 0 0 1;0 -1 -2 4;0 -2 1 3;1 4 3 1];hohd(a);

得到原始实数矩阵的拟上三角矩阵A。

A =

2.0000 -1.0000 0.0000 0.0000

-1.0000 1.0000 5.0000 -0.0000

0.0000 5.0000 -2.2000 -0.4000

0.0000 -0.0000 -0.4000 2.2000

2、执行命令：

>> clear,clc,A=[2 -1 0 0;-1 1 5 0; 0 5 -2.2 -0.4; 0 0 -0.4 2.2];qrfenjie(A);

得到迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。

i =

38

Ak =

-5.9068 0.0315 -0.0000 0.0000

0.0315 4.8972 -0.0000 -0.0000

0.0000 -0.0000 2.2138 0.0007

0.0000 0.0000 0.0007 1.7958

lanmda =

-5.9068 4.8972 2.2138 1.7958

3、执行命令：

>> clear,clc,A=[2 -1 0 0;-1 1 5 0; 0 5 -2.2 -0.4; 0 0 -0.4 2.2];dp_qrfenjie(A);

得到用带位移的QR分解法所用的迭代次数i，舒尔矩阵Ak，以及A的特征值矩阵lamda。

i =

8

Ak =

-5.9068 -0.0056 -0.0000 -0.0000

-0.0056 4.8973 0.0000 0.0000

-0.0000 0.0000 1.7958 0.0000

0.0000 0.0000 0.0000 2.2138

lanmda =

-5.9068 4.8973 1.7958 2.2138

4、执行命令：

>> clear,clc,a=[2 0 0 1;0 -1 -2 4;0 -2 1 3;1 4 3 1],p=[-5.9068 4.8972 2.2138 1.7958];vect(a,p);

得到Z矩阵，Z矩阵每一行的第一个数是迭代次数，第二个数是特征值，后面的数为特征向量。

Z =

Columns 1 through 4

5.000000000000000 -5.906847942119164 0.112419521552950 1.000000000000000

5.000000000000000 4.897302326740487 0.345148654584838 0.505132030845888

5.000000000000000 2.213757601733807 -0.282974649968150 -0.697610774659887

5.000000000000000 1.795788013644870 1.000000000000000 -0.324049********* Columns 5 through 6

0.675655845769651 -0.888884062644966

0.510541849590699 1.000000000000000

1.000000000000000 -0.060487982528655

0.044562197436887 -0.204211986355130

六、结果分析：

1、QR算法是求取一般矩阵所有特征值的一种非常有效的数值方法。

2、由2、3两步的运行结果看出带位移的QR算法与一般的QR算法相比具有更好的收敛效

果，极大地减少了运算量。

3、通过反幂法只需比较小的计算量就可以得到精度更高的特征值，以及对应的特征向量。

数值分析大作业-三、四、五、六、七

大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0; end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +??=-= ?-???解：Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while nerrorlim n=n+1; else break ; end x0=x; end disp(['迭代次数: n=',num2str(n)]) disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) （2）子函数 非线性函数f function y=f(x) y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); end

数值计算方法比较

有限差分方法(FDM:Finite Difference Method)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。有限差分法主要集中在依赖于时间的问题（双曲型和抛物型方程）。有限差分法方面的经典文献有Richtmeyer & Morton的《Difference Methods for Initial-Value Problems》；R. LeVeque《Finite Difference Method for Differential Equations》；《Numerical Methods for C onservation Laws》。 注：差分格式： （1）从格式的精度来划分，有一阶格式、二阶格式和高阶格式。 （2）从差分的空间形式来考虑，可分为中心格式和逆风格式。 （3）考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。 目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。 构造差分的方法： 构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。 有限差分法的不足：由于采用的是直交网格，因此较难适应区域形状的任意性，而且区分不出场函数在区域中的轻重缓急之差异，缺乏统一有效的处理自然边值条件和内边值条件的方法，难以构造高精度(指收敛阶)差分格式，除非允许差分方程联系更多的节点(这又进一步增加处理边值条件韵困难)。另外它还有编制不出通用程序的困难。 有限差分法的优点：该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念 直观，表达简单，精度可选而且在一个时间步内，对于一个给定点来说其相关的空间点只是 与该相邻的几点，而不是全部的空间点。是发展较早且比较成熟的数值方法 广义差分法（有限体积法）（GDM:Generalized Difference Method）：1953年，Mac—Neal 利用积分插值法(也称积分均衡法)建立了三角网格上的差分格 式，这就是以后通称的不规划网格上的差分法．这种方法的几何误差小，特别是给出了处理自然边值条件(及内边值条件)的有效方法，堪称差分法的一大进步。1978年，李荣华利用有限元空间和对偶单元上特征函数的推广——局部Taylor展式的公项，将积分插值法改写成广义Galerkin法形式，从而将不规则网格差分法推广为广义差分法.其基本思路是，将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有

数值分析大作业三 四 五 六 七

大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:');

flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0;

end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =-= ?-???解： Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1;

数值分析(计算方法)总结

第一章绪论 误差来源：模型误差、观测误差、截断误差（方法误差）、舍入误差 是的绝对误差，是的误差，为的绝对误差限（或误差限） 为的相对误差，当较小时，令 相对误差绝对值得上限称为相对误差限记为：即： 绝对误差有量纲，而相对误差无量纲 若近似值的绝对误差限为某一位上的半个单位，且该位直到的第一位非零数字共有n位，则称近似值有n位有效数字，或说精确到该位。 例：设x==…那么,则有效数字为1位，即个位上的3，或说精确到个位。 科学计数法：记有n位有效数字，精确到。 由有效数字求相对误差限：设近似值有n位有效数字，则其相对误差限为 由相对误差限求有效数字：设近似值的相对误差限为为则它有n位有效数字 令 1.x+y近似值为和的误差（限）等于误差（限）的 和 2.x-y近似值为 3.xy近似值为 4. 1．避免两相近数相减 2．避免用绝对值很小的数作除数 3．避免大数吃小数 4．尽量减少计算工作量 第二章非线性方程求根 1.逐步搜索法 设f (a) <0, f (b)> 0，有根区间为(a, b)，从x0=a出发，按某个预定步长(例如h=(b-a)/N)

一步一步向右跨，每跨一步进行一次根的搜索，即判别f(x k)=f(a+kh)的符号，若f(x k)>0(而 f(x k-1)<0),则有根区间缩小为[x k-1,x k] (若f(x k)=0，x k即为所求根), 然后从x k-1出发，把搜索步长再缩小，重复上面步骤，直到满足精度：|x k-x k-1|0.将[a0,b0]对分，中点x0= ((a0+b0)/2),计算 f(x0)。 3.比例法 一般地，设[a k,b k]为有根区间，过(a k, f(a k))、(b k, f(b k))作直线，与x轴交于一点x k,则： 1.试位法每次迭代比二分法多算一次乘法，而且不保证收敛。 2.比例法不是通过使求根区间缩小到0来求根，而是在一定条件下直接构造出一个点列（递推公式），使该点列收敛到方程的根。——这正是迭代法的基本思想。 事先估计: 事后估计 局部收敛性判定定理： 局部收敛性定理对迭代函数的要求较弱，但对初始点要求较高，即初始点必须选在精确解的附近 Steffensen迭代格式： Newton法： Newton下山法：是下山因子 弦割法： 抛物线法：令 其中：

数值分析作业答案

数值分析作业答案 插值法 1、当x=1，-1，2时，f(x)=0，-3，4，求f(x)的二次插值多项式。 （1）用单项式基底。 （2）用Lagrange插值基底。 （3）用Newton基底。 证明三种方法得到的多项式是相同的。 解：（1）用单项式基底 设多项式为: ， 所以： 所以f(x)的二次插值多项式为： （2）用Lagrange插值基底 Lagrange插值多项式为： 所以f(x)的二次插值多项式为： (3) 用Newton基底: 均差表如下： xk f(xk) 一阶均差二阶均差 1 0 -1 -3 3/2 2 4 7/ 3 5/6 Newton插值多项式为： 所以f(x)的二次插值多项式为： 由以上计算可知，三种方法得到的多项式是相同的。 6、在上给出的等距节点函数表，若用二次插值求ex的近似值，要使截断误差不超过10-6，问使用函数表的步长h应取多少？ 解：以xi-1,xi,xi+1为插值节点多项式的截断误差，则有 式中 令得 插值点个数

是奇数，故实际可采用的函数值表步长 8、，求及。 解：由均差的性质可知，均差与导数有如下关系： 所以有： 15、证明两点三次Hermite插值余项是 并由此求出分段三次Hermite插值的误差限。 证明：利用[xk，xk+1]上两点三次Hermite插值条件 知有二重零点xk和k+1。设 确定函数k(x): 当或xk+1时k(x)取任何有限值均可； 当时，，构造关于变量t的函数 显然有 在[xk,x][x,xk+1]上对g(x)使用Rolle定理，存在及使得 在，，上对使用Rolle定理，存在，和使得 再依次对和使用Rolle定理，知至少存在使得 而，将代入，得到 推导过程表明依赖于及x 综合以上过程有： 确定误差限： 记为f(x)在[a,b]上基于等距节点的分段三次Hermite插值函数。在区间[xk，xk+1]上有 而最值 进而得误差估计： 16、求一个次数不高于4次的多项式，使它满足，，。

数值计算方法大作业

目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20)

4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28)

数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10）

1.2牛顿法 程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10）

数值分析大作业

数值分析报大作业 班级：铁道2班 专业：道路与铁道工程 姓名：蔡敦锦 学号：13011260

一、序言 该数值分析大作业是通过C语言程序编程在Microsoft Visual C++ 6.0编程软件上运行实现的。本来是打算用Matlab软间来计算非线性方程的根的。学习Matlab也差不多有一个多月了，感觉自己编程做题应该没什么问题了；但是当自己真心的去编程、运行时才发现有很多错误，花了一天时间修改、调试程序都没能得到自己满意的结果。所以，我选择了自己比较熟悉的C程序语言来编程解决非线性的求值问题，由于本作业是为了比较几种方法求值问题的收敛速度和精度的差异，选择了一个相对常见的非线性函数来反映其差异，程序运行所得结果我个人比较满意。编写C语言，感觉比较上手，程序出现问题也能比较熟练的解决。最终就决定上交一份C程序语言编程的求值程序了！

二、选题 本作业的目的是为了加深对非线性方程求根方法的二分法、简单迭代法、、牛顿迭代法弦截法等的构造过程的理解；能将各种方法的算法描述正确并且能够改编为程序并在计算机上实现程序的正确合理的运行，能得到自己满意的结果，并且能调试修改程序中可能出现的问题和程序功能的增减修改。本次程序是为了比较各种方法在求解同一非线性方程根时，在收敛情况上的差异。 为了达到上面的条件我选择自己比较熟悉的语言—C语言来编程，所选题目为计算方程f(x)=x3-2x-5=0在区间[2,3]内其最后两近似值的差的绝对值小于等于5 ?的根的几种方法的比较。 110- 本文将二分法、牛顿法、简单迭代法、弦截法及加速收敛法这五种方法在同一个程序中以函数调用的方式来实现，比较简洁明了，所得结果能很好的比较，便于分析；发现问题和得出结论。

数值分析第一次作业及参考答案

数值计算方法第一次作业及参考答案 1. 已测得函数()y f x =的三对数据：（0，1），（－1，5），（2，－1）， （1）用Lagrange 插值求二次插值多项式。（2）构造差商表。（3）用Newton 插值求二次插值多项式。 解：(1)Lagrange 插值基函数为 0(1)(2)1 ()(1)(2)(01)(02)2 x x l x x x +-= =-+-+- 同理 1211 ()(2),()(1)36 l x x x l x x x = -=+ 故 2 20 2151 ()()(1)(2)(2)(1) 23631 i i i p x y l x x x x x x x x x =-==-+-+-++=-+∑ （2）令0120,1,2x x x ==-=，则一阶差商、二阶差商为 011215 5(1) [,]4, [,]20(1) 12 f x x f x x ---= =-= =----- 0124(2) [,,]102 f x x x ---= =- 实际演算中可列一张差商表： （3）用对角线上的数据写出插值多项式 2 2()1(4)(0)1*(0)(1)31P x x x x x x =+--+-+=-+ 2. 在44x -≤≤上给出()x f x e =的等距节点函数表，若用二次插值求x e 的近似值，要使 截断误差不超过6 10-，问使用函数表的步长h 应取多少 解： ()40000(), (),[4,4],,,, 1.x k x f x e f x e e x x h x x h x x th t ==≤∈--+=+≤考察点及

(3) 2000 4 43 4 3 () ()[(()]()[()] 3! (1)(1) (1)(1) 3!3! .(4,4). 6 f R x x x h x x x x h t t t e t h th t h e h e ξ ξ =----+ -+ ≤+??-= ≤∈- 则 4 36 ((1)(1) 100.006. t t t h - -+± << Q在点 得 3.求2 () f x x =在［a,b］上的分段线性插值函数() h I x，并估计误差。 解： 22 22 11 1 111 22 11 11 1 () () k k k k h k k k k k k k k k k k k k k k k k k x x x x x x I x x x x x x x x x x x x x x x x x x x x x ++ + +++ ++ ++ + --- =+= --- ?-? -=+- - [] 2 11 22 11 ()()()[()] 11 ()() 44 h h k k k k k k k k R x f x I x x x x x x x x x x x x x h ++ ++ =-=-+- =--≤-= 4.已知单调连续函数() y f x =的如下数据 用插值法计算x约为多少时() 1. f x=（小数点后至少保留4位） 解：作辅助函数()()1, g x f x =-则问题转化为x为多少时，()0. g x=此时可作新 的关于() i g x的函数表。由() f x单调连续知() g x也单调连续，因此可对() g x的数值进行反插。的牛顿型插值多项式为 1()0.110.097345( 2.23)0.451565( 2.23)( 1.10) 0.255894( 2.23)( 1.10)(0.17) x g y y y y y y y - ==-+++++ -++-

数值分析作业答案part

6.4.设??? ? ? ??=5010010a b b a A ，0det ≠A ，用a ，b 表示解线性方程组f Ax =的雅可比迭代与 高斯—塞德尔迭代收敛的充分必要条件。 解 雅可比迭代法的迭代矩阵 ? ??? ??? ? ??----=???? ? ??----????? ??=-050100100100000001010101 a b b a a b b a B J ， ?? ? ?? -=-1003||2ab B I J λλλ，10||3)(ab B J = ρ。 雅可比迭代法收敛的充分必要条件是3 100 ||

数值计算方法》试题集及答案

《计算方法》期中复习试题 一、填空题： 1、已知3.1)3(,2.1)2(,0.1)1(===f f f ，则用辛普生（辛卜生）公式计算求得 ?≈3 1 _________ )(dx x f ,用三点式求得≈')1(f 。 答案：2.367，0.25 2、1)3(,2)2(,1)1(==-=f f f ，则过这三点的二次插值多项式中2 x 的系数为 ，拉 格朗日插值多项式为 。 答案：-1， )2)(1(21 )3)(1(2)3)(2(21)(2--------= x x x x x x x L 3、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字； 4、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( )； 答案 )(1)(1n n n n n x f x f x x x '--- =+ 5、对1)(3 ++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 )； 6、计算方法主要研究( 截断 )误差和( 舍入 )误差； 7、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时，二分n 次后的误差限为 ( 1 2+-n a b )； 8、已知f (1)＝2，f (2)＝3，f (4)＝5.9，则二次Newton 插值多项式中x 2系数为( 0.15 )； 11、 两点式高斯型求积公式?1 d )(x x f ≈( ?++-≈1 )] 321 3()3213([21d )(f f x x f )，代数精度 为( 5 )； 12、 为了使计算 32)1(6 )1(41310-- -+-+ =x x x y 的乘除法次数尽量地少，应将该表达 式改写为 11 ,))64(3(10-= -++=x t t t t y ，为了减少舍入误差，应将表达式1999 2001-

北航数值分析报告第三次大作业

数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： x y当作已知量代入题目给定的非线性方程组，求（1）解非线性方程组。将给定的(,) i i

得与(,)i i x y 相对应的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和(,)i i p x y 的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解 (,)i i p x y 可以直接使用（3）中的C[r][s]和k 。 （二）具体算法设计： （1）解非线性方程组 牛顿法解方程组()0F x =的解* x ，可采用如下算法： 1）在* x 附近选取(0) x D ∈，给定精度水平0ε>和最大迭代次数M 。 2）对于0,1, k M =执行 ① 计算() ()k F x 和()()k F x '。 ② 求解关于() k x ?的线性方程组 () ()()()()k k k F x x F x '?=- ③ 若() () k k x x ε∞∞ ?≤，则取*()k x x ≈，并停止计算；否则转④。 ④ 计算(1) ()()k k k x x x +=+?。 ⑤ 若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。 （2）分片双二次插值 给定已知数表以及需要插值的节点，进行分片二次插值的算法： 设已知数表中的点为： 00(0,1,,) (0,1,,)i j x x ih i n y y j j m τ=+=???=+=?? ，需要插值的节点为(,)x y 。 1) 根据(,)x y 选择插值节点(,)i j x y ： 若12h x x ≤+ 或12 n h x x ->-，插值节点对应取1i =或1i n =-，

北航数值分析大作业第二题精解

目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。

这样即求出所有A所有实特征值对应的一个特征向量。 #include #include #include #define N 10 #define E 1.0e-12 #define MAX 10000 //以下是符号函数 double sgn(double a) { double z; if(a>E) z=1; else z=-1; return z; } //以下是矩阵的拟三角分解 void nishangsanjiaodiv(double A[N][N]) { int i,j,k; int m=0; double d,c,h,t; double u[N],p[N],q[N],w[N]; for(i=0;i

上海大学_王培康_数值分析大作业

数值分析大作业（2013年5月） 金洋洋（12721512），机自系 1.下列各数都是经过四舍五入得到的近似值，试分别指出它 们的绝对误差限, 相对误差限和有效数字的位数。 X1 =5.420, x 2 =0.5420, x 3=0.00542, x 4 =6000, x 5=50.610? 解：根据定义：如果*x 的绝对误差限 不超过x 的某个数位的半个单位，则从*x 的首位非零数字到该位都是有效数字。 显然根据四舍五入原则得到的近视值，全部都是有效数字。 因而在这里有：n1=4, n2=4, n3=3, n4=4, n5=1 （n 表示x 有效数字的位数） 对x1：有a1=5, m1=1 （其中a1表示x 的首位非零数字，m1表示x1的整数位数） 所以有绝对误差限 143 11 (1)101022 x ε--≤ ?=? 相对误差限 31() 0.510(1)0.00923%5.4201 r x x x εε-?= == 对x2：有a2=5, m2=0 所以有绝对误差限 044 11 (2)101022 x ε--≤ ?=? 相对误差限 42() 0.510(2)0.00923%0.54202 r x x x εε-?= == 对x3：有a3=5, m3=-2 所以有绝对误差限 235 11 (3)101022 x ε---≤ ?=? 相对误差限 53() 0.510(3)0.0923%0.005423 r x x x εε-?= == 对x4：有a4=0, m4=4 所以有绝对误差限 4411(4)1022 x ε-≤?= 相对误差限 4() 0.5 (4)0.0083%6000 4 r x x x εε= = = 对x5：有a5=6, m5=5 所以有绝对误差限 514 11(5)101022 x ε-≤ ?=? 相对误差限 45() 0.510(5)8.3%600005 r x x x εε?= ==

数值分析大作业QR分解

题目:设计求取n n ?实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵 20010-1-24A=0-2131 43 1?? ? ? ? ??? 为例进行具体的求解计算。 一、 算法分析： 一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。 二、 算法设计： 1、 原始实矩阵A R n n ?∈拟上三角化 为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’ 记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)r ij a i n j r r n ==+。 对于r=1,2,…,n -2执行 (1) 若() (2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。 (2) 计算 令 ()2 ()() 1,1,1,sgn(0,sgn()=1) r r r r r r r r r r r c a a a ρ+++=-=，（若则取 (3) 令-0=r n r u u ?? ? ?? ，()()()-1,2,1(,,...,)r r r T n r r r r r r nr u a c a a ρ++=- (4) 计算 (r)(r)(r)T n-r r+1,r r r+2,r nr r T n-r T n-r n-r n-r n-r r+1r 1u = (a -c ,a ,...,a )ρ I H =I -2uu =H H =I -2u u A =HA H ?? ? ?? (5) 继续 算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。 2、 拟上三角矩阵QR 分解的求原矩阵的全部特征值 记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵，

数值分析作业答案(第5章)

5.1.设A 是对称矩阵且011≠a ，经过一步高斯消去法后，A 约化为 ?? ????21 110 A a a T 证明2A 是对称矩阵。 证明 由消元公式及A 的对称性，有 ,,,3,2,,)2(111 11111 )2(n j i a a a a a a a a a a ji i j ji j i ij ij ==-=- = 故2A 对称。 5.2.设n ij a A )(=是对称正定矩阵，经过高斯消去法一步后，A 约化为 ?? ????21 110 A a a T 其中1)2(2)(-=n ij a A 。证明： (1).A 的对角元素;,,2,1,0n i a ii => (2).2A 是对称正定矩阵。 证明 (1).因为A 对称正定，所以 n i e Ae a i i ii ,,2,1,0),( =>=， 其中T i e )0,,0,1,0,,0( =为第i 个单位向量。 (2).由A 的对称性及消元公式，有 ,,,3,2,,)2(111 11111 )2(n j i a a a a a a a a a a ji i j ji j i ij ij ==-=- = 故2A 也对称。 又由A L A a a T 121110=????? ?，其中

??? ?????- =? ????? ? ?????????--=-111 1 11111 21101 1011n n I a a a a a a L ， 可见1L 非奇异，因而对任意0≠x ，由A 的正定性，有 ,0),(),(,011111>=≠x AL x L x AL L x x L T T T T 故T AL L 11正定。 由,000110211 111121111 1?? ? ?? ?=????????-??????=-A a I a a A a a AL L n T T T 而011>a ，故知2A 正定

数值分析大作业 三、四、五、六、七

大作业 三 1. 给定初值0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0; end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =- = ?-???解：Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while nerrorlim n=n+1; else break ; end x0=x; end disp(['迭代次数: n=',num2str(n)]) disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) （2）子函数 非线性函数f function y=f(x) y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); end

数值分析作业答案.doc

第2章 插值法 1、当x=1，-1，2时，f(x)=0，-3，4，求f(x)的二次插值多项式。 （1）用单项式基底。 （2）用Lagrange 插值基底。 （3）用Newton 基底。 证明三种方法得到的多项式是相同的。 解：（1）用单项式基底 设多项式为:2 210)(x a x a a x P ++=， 所以：64 211111 1111122 2 211 200 -=-==x x x x x x A 3 76144 211111114241 13110111)() ()(22 221120 022 2 22 11 120 00-=-= ---==x x x x x x x x x f x x x f x x x f a 2 3694211111114411 31101111)(1)(1 )(122 221120 02 2 22112 001=--= --==x x x x x x x x f x x f x x f a 6 5654 2 1 1111114 2 1 3 11011111) (1)(1)(122 2 21120 022 11 00 2=--= ---==x x x x x x x f x x f x x f x a 所以f(x)的二次插值多项式为：26 52337)(x x x P ++-= （2）用Lagrange 插值基底 )21)(11() 2)(1())(())(()(2010210-+-+=----=x x x x x x x x x x x l )21)(11() 2)(1())(())(()(2101201------=----=x x x x x x x x x x x l ) 12)(12() 1)(1())(())(()(1202102+-+-=----= x x x x x x x x x x x l

北航数值分析报告大作业第八题

北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期2014 年11月26 日

一．题目 关于x , y , t , u , v , w 的方程组(A.3) ???? ?? ?=-+++=-+++=-+++=-+++79 .0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。 表A-1 二维数表 t z u 0 0.4 0.8 1.2 1.6 2 0 -0.5 -0.34 0.14 0.94 2.06 3.5 0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38 0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42 0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62 0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02 1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5 1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式 ∑∑===k i k j s r rs y x c y x p 00 ),( 要求p (x , y )以最小的k 值达到以下的精度 ∑∑==-≤-=10020 7210)],(),([i j i i i i y x p y x f σ 其中j y i x i i 05.05.0,08.0+==。 2. 计算),(),,(* ***j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼 近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。

数值分析计算方法试题集及答案

数值分析复习试题 第一章 绪论 一. 填空题 1.* x 为精确值 x 的近似值；() **x f y =为一元函数 ()x f y =1的近似值； ()**,*y x f y =为二元函数()y x f y ,2=的近似值，请写出下面的公式：**e x x =-： *** r x x e x -= ()()()*'1**y f x x εε≈? ()() () ()'***1**r r x f x y x f x εε≈ ? ()()()() ()* *,**,*2**f x y f x y y x y x y εεε??≈?+??? ()()()()() ** * *,***,**222r f x y e x f x y e y y x y y y ε??≈ ?+??? 2、 计算方法实际计算时，对数据只能取有限位表示，这时所产生的误差叫 舍入误 差 。 3、 分别用2.718281，2.718282作数e 的近似值，则其有效数字分别有 6 位和 7 位；又取 1.73≈-21 1.73 10 2 ≤?。 4、 设121.216, 3.654x x ==均具有3位有效数字，则12x x 的相对误差限为 0.0055 。 5、 设121.216, 3.654x x ==均具有3位有效数字，则12x x +的误差限为 0.01 。 6、 已知近似值 2.4560A x =是由真值T x 经四舍五入得 到,则相对误差限为 0.0000204 . 7、 递推公式,??? ? ?0n n-1y =y =10y -1,n =1,2, 如果取0 1.41y ≈作计算,则计算到10y 时,误 差为 81 10 2 ?;这个计算公式数值稳定不稳定 不稳定 . 8、 精确值 14159265.3* =π，则近似值141.3*1=π和1415.3*2=π分别有 3