北科大研究生计算方法作业

计算方法

姓名：

学号：

班级：

指导教师：

目录

作业1 (1)

作业2 (5)

作业3 (8)

作业4 (10)

作业5 (14)

作业6 (16)

作业7 (17)

作业1

1、分别用不动点迭代与Newton 法求解方程 -+=x 2x e 30的正根与负根。 解：

(1)不动点迭代

a.原理：

将 230x x e -+=变型为1()k k x g x +=进行迭代，直到为止 变型后为有两种形式：和 b.程序：初值为1

形式:

x=zeros(100,1); tol=1; i=1; x(1)=1;

whiletol>=10e-6; disp(x(i))

x(i+1)=log(2*x(i)+3); tol=abs(x(i+1)-x(i)); i=i+1; end

disp(i-1); 形式:

x=zeros(100,1); tol=1; i=1; x(1)=1;

whiletol>=10e-6; disp(x(i))

x(i+1)=(exp(x(i))-3)/2; tol=abs(x(i+1)-x(i)); i=i+1; end disp(i-1);

c.运行结果：

初值为1

(23)

1ln

k x k x ++=6

110k k x x -+-<13

2

k x k e x +-=

(23)1ln k x k x ++=132

k x k e x +-=

迭代次数：11

迭代次数：9

(2)Nexton法

a.原理：

令

()

()

1'

k

k k

k

f x

x x

f x

+

=-得到迭代公式为：

()

1

23

2

k

k

x

k

k k x

x e

x x

e

+

-+

=-

-

b.程序：初值为0

x=zeros(100,1);

tol=1;

i=1;

x(1)=0;

whiletol>=10e-6;

disp(x(i))

x(i+1)=x(i)-((2*x(i)-exp(x(i))+3)/(2-exp(x(i))));

tol=abs(x(i+1)-x(i));

i=i+1;

end

disp(i-1);

初值为1

x=zeros(100,1);

tol=1;

i=1;

x(1)=1;

whiletol>=10e-6;

disp(x(i))

x(i+1)=x(i)-((2*x(i)-exp(x(i))+3)/(2-exp(x(i))));

tol=abs(x(i+1)-x(i));

i=i+1;

end

disp(i-1)

a=x(i-1);

b=2*a-exp(a)+3;

disp(b);

c.运行结果：

初值为0

迭代次数：5

初值为1

迭代次数：8 -1.6171e-006

结果分析：

不动点迭代会因为迭代公式选取的不同得出不同的迭代结果，而牛顿法迭代会因为初值选取的不同而得到不同的结果。牛顿法比不动点迭代法收敛速度快，能较少的迭代达到理想的结果。

2、用Newton 法求解方程-x sinx 0＝的根.再用Steffensen’smethod 加速其收敛。 解： (1)Newton 法 a.原理：

令()()

1'

k k k k f x x x f x +=-

得到迭代公式为：1sin 1cos k k

k k k x x x x x +-=-- b.程序：

x=zeros(100,1); tol=1; i=1; x(1)=1;

whiletol>=10e-4; disp(x(i))

x(i+1)=x(i)-((x(i)-sin(x(i)))/(1-cos(x(i)))); tol=abs(x(i+1)-x(i)); i=i+1; end

disp(i-1) a=x(i-1); b=a-sin(a); disp(b);

迭代次数：16

1.8120e-009

(2)Steffensen’smethod加速

a.原理：

Steffensen加速是Aitken加速与不动点迭代的结合

b.程序：

p0=5;N=20;tol=10e-6;

n=0;p(1)=p0;

while n<=N

for k=1:2

p(k+1)=p(k)-((p(k)-sin(p(k)))/(1-cos(p(k))));

end

p1=p(1)-(p(2)-p(1))^2/(p(3)-2*p(2)+p(1));

f0=p1-sin(p1);

if abs(f0)

break

end

n=n+1;p(1)=p1;

end

disp(p1);disp(n)

c.运行结果：

-3.213248067584773e-004

迭代次数：2

结果分析：

Steffensen加速极大的改善了原迭代的收敛性质。用牛顿迭代16次，使用Steffensen加速后，仅需迭代2次即得到满足要求的迭代结果。

作业2

1、分别用Jacobi 、Seidel 、Sor(w=1.1,1.2,1.3,1.4,1.5)方法求解方程组x =A b ，

这里1010121111211112A ?-????-??=????-?? ，101

11b ?-??

??

??

=????

-?? []000T

X = ；510e -= 解：

(1)Jacobi 迭代

a.原理：

Jacobi 迭代方程：

()11()()11,

11,2,,0,1,i n

k k k i i ij j ij j j j i x b a x a x i n k aii -+==+??=--== ???∑∑

b.程序：

A=ones(10); fori=1:10; A(i,i)=-12; end

b=ones(10,1); B=-1*b;

X0=zeros(10,1); Tol=10^-6; N=1000; X=X0;

for K=1:N fori=1:10

X(i)=(B(i)-A(i,:)*X0)/A(i,i)+X0(i); if norm(X-X0)

X0=X; end

disp('发散') c.运行结果：

迭代次数：41

(2)Seidel 迭代 a.原理：

seidel 迭代方程：

()11(1)()11,

11,2,,0,1,i n

k k k i i ij j ij j j j i x b a x a x i n k aii -++==+??=--== ???∑∑ b.程序：

A=ones(10); fori=1:10; A(i,i)=-12; end

b=ones(10,1); B=-1*b;

X0=zeros(10,1); Tol=10^-6; N=1000; X=X0;

for K=1:N; fori=1:10

X(i)=(B(i)-A(i,:)*X)/A(i,i)+X(i); if norm(X-X0)

disp('发散') c.运行结果：

迭代次数：24

(3)Sor 迭代 a.原理：

sor 迭代方程：

()

11()(1)()11,

(1)1,2,,0,1,i n

k k k k i

i

i ij j ij i j j i w x w x

b a x a x i n k aii -++==+??=-+--== ???∑∑ b.程序：

A=ones(10); fori=1:10; A(i,i)=-12; end

b=ones(10,1);

B=-1*b;

X0=zeros(10,1);

w=1.1;

Tol=10^-6;

N=1000;

X=X0;

for K=1:N;

fori=1:10

X(i)=w*(B(i)-A(i,:)*X)/A(i,i)+X(i);

if norm(X-X0)

disp(X);disp(K);

return

end

end

X0=X;

end

disp('发散')

c.运行结果：

迭代次数：19

结果分析：

对比三种迭代方法的结果和迭代次数，初步得出结论：若Jacobi迭代法和

Gauss-Seidel迭代法都收敛时，Gauss-Seidel迭代法比用于Jacobi迭代法求解收敛速度快。应用Sor时,选取合适的松弛因子是很关键的，从解该方程组程序及结果对照表中可以看出，松弛因子w的选取，可以影响其收敛速度。

作业3

1、用Newton 法与最速下降法求方程组

222

220x y x y

+-==

在(0.8,0.7)附近的根。 解：

(1)Newton 法

a.原理：

由题目所给方程组可以得到

2222224(),'()21x y x y F x F x x x y ??+-??== ? ? ?--????

取初始向量()0(0.8,0.7)T x =，应用牛顿迭代法

()()()()()()1(1)()'(,)(,)0,1,2k k k k k k k k F x y z F x y k x x z y y ++?=-?

=?????

=+ ? ??

?????

b.程序：

x0=[0.8;0.7]; Tol=10^-6; N=1000; for K=1:N

F=[x0(1)^2+2*x0(2)^2-2;x0(1)^2-x0(2)]; dF=[2*x0(1),4*x0(2);2*x0(1),-1]; y=-dF\F; x=x0+y;

if norm(y)

disp('发散') c.运行结果：

x=[0.8836;0.7808] 迭代次数：4

(2)最速下降法

a.原理：

由原方程组构造一个模函数 ()2212()()x f x f x φ=+

其中

()()2221222,f x x y f x x y =+-=-

b.程序：

function h=h(t) X=[0.8; 0.7];

dF=[2*(X(1)^2+2*X(2)^2-2)*2*X(1)+2*(X(1)^2-X(2))*2*X(1

);2*(X(1)^2+2*X(2)^2-2)*4*X(2)+2*(X(1)^2-X(2))*(-1)];

X1=X-t*dF;

h=(X1(1)^2+2*X1(2)^2-2)^2+(X1(1)^2-X1(2))^2;

X0=[0.8;0.7];Tol=10^-6;N=1000; for K=1:N

F=(X0(1)^2+2*X0(2)^2-2)^2+(X0(1)^2-X0(2))^2;

dF=[2*(X0(1)^2+2*X0(2)^2-2)*2*X0(1)+2*(X0(1)^2-X0(

2))*2*X0(1);2*(X0(1)^2+2*X0(2)^2-2)*4*X0(2)+2*(X0(1)^2-X0(2))*(-1)];

U=fminbnd('h',0,100); X0=X0-U*dF; if F

disp(X0);disp(K); return; end end c.运行结果：

x=[0.8834;0.7809] 迭代次数：11

结果分析：

最速下降法计算简便，但其收敛速度是线性收敛，不如牛顿法快，它的优点是对任意初值(0)x 都收敛，最速下降法非常可靠。

作业4

1、用幂法与反幂法求矩阵A 的按模最大、最小特征值与对应的特征向量。

411113*********

2A ????-?

?=??-?

???

解：

(1)幂法： a.原理：

求解最大特征值与特征向量

b.程序：

A=[4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2]; v=[1;1;1;1]; r=10; err=1; k=1; N=100;

while(k<=N && err>eps) u=A*v;

[m j]=max(abs(u)); dc=abs(r-m); u=u/m;

dv=norm(u-v); err=max(dc,dv); v=u; r=m; k=k+1; end r u c.运行结果：

r =5.2361 u =1.0000 0.6180 0.1180

0.5000

(2)反幂法：

a.原理：

求解最小特征值和特征向量 反幂法的计算公式：

()()()()

()()()

11

/max 0,1,2,k k k k k U V V k V

A qI U -+?=?

=??=-? b.程序：

A=[4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2]; A1=inv(A); v=[1;1;1;1]; r=10;

err=1;k=1;N=100;

while(k<=N && err>eps) u=A1*v;

[m j]=max(abs(u)); dc=abs(r-m); u=u/m;

dv=norm(u-v); err=max(dc,dv); v=u; r=m; k=k+1; end r=1/r; r u

c.运行结果：

r =0.7639 u =-0.4721 0.7639 1.0000 -0.2361

2、用House-holder 变换求矩阵A 的QR 分解，并用QR 方法做三次迭代。

4111131111201102A ????

-?

?=??-????

解：

(1)QR 分解 a.原理：

Householder 变换公式：

2T H I uu =-

利用Householder 矩阵求矩阵A 的QR 分解公式：

1231121****0*...0*0110n n n c c c A H H H A c --????

??

??==????????

121n Q H H H -= ，1n R A -=

b.程序：

A=[4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2]; n=size(A); Q=eye(n); for k=1:n-1

c=norm(A(k:n,k)); w=zeros(n,1);

w(k+1:n)=A(k+1:n,k);

w(k)=A(k,k)+sign(A(k,k))*c; P=eye(n)-2*w*w'/norm(w)^2; Q=Q*P; A=P*A; end R=A;

disp(Q); disp(R); c.运行结果：

0.91770.1543

0.33620.14510.22940.85740.35790.2902 0.22940.44580.86220.07250.22940.20580.12470.9431Q --????----??

=?

?

---????---??4.35891.60591.14711.60590.0000 3.06941.90341.11460.00000.00001.03030.27110.0000

0.00000.0000

1.4510R ----????--??

=?

?

--??

????

(2)QR 方法做三次迭代 a.原理： b.程序：

A=[4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2]; n=size(A); Q=eye(n);

N=1;

while N<=3 for k=1:n-1

c=norm(A(k:n,k)); w=zeros(n,1); w(k+1:n)=A(k+1:n,k);

w(k)=A(k,k)+sign(A(k,k))*c;

P=eye(n)-2*w*w'/norm(w)^2; Q=Q*P; A=P*A; end

A=A*Q; N=N+1; end

R=A*inv(Q); disp(Q); disp(R);

c.运行结果：

0.80370.28340.41130.32340.50460.63420.58410.04480.13220.71930.67450.10100.28640.01000.18620.9398Q --????----??=??

--????--?? 4.68410.03231.62731.58020.0000 3.15971.05091.25590.00000.00000.85060.015600.00000.00001.5887R ---????--??=??--??????

作业5

1、目的：观察lagrange 插值的Runge 现象。 对于函数2

1

()125x

f x =

+进行lagrange 插值。取不同的等分数n(n=5,n=10)，将区间[-1，1]n 等分，取等距节点。把插值多项式的曲线画在同一张图上进行比较。 a.原理：

朗格朗日插值 b.程序：原曲线

function a=f1(b) a=1/(1+25*b^2);

n=5插值曲线

function y1=f2(x1) n=5; h=2/n; x(1)=-1;

y(1)=f1(x(1)); fori=1:n

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

y(i+1)=f1(x(i+1)); end Y=0;

for k=1:n+1 t=1; fori=1:n+1

ifi~=k; t=t*(x1-x(i))/(x(k)-x(i)); end end

Y=Y+t*y(k); end y1=Y;

n=10插值曲线

function y1=f3(x1) n=10; h=2/n; x(1)=-1;

y(1)=f1(x(1)); fori=1:n

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

y(i+1)=f1(x(i+1)); end

Y=0;

for k=1:n+1

t=1;

fori=1:n+1

ifi~=k;t=t*(x1-x(i))/(x(k)-x(i));

end

end

Y=Y+t*y(k);

end

y1=Y;

n=1000;

h=2/n;

fori=1:1000

a=-1+(i-1)*h;

t(i)=-1+(i-1)*h;

y1(i)=f2(a);

y2(i)=f3(a);

y3(i)=f1(a);

end

plot(t,y1,'r-..',t,y2,'b-.',t,y3,'g-')

c.运行结果：

结果分析：

利用lagrange插值拟合曲线并不是插值点越多越接近真实曲线。

1、用Romberg 求积公式计算积分：3

1

?x e cosxdx ,精度达到610-。

a.原理：

Romberg 求积公式：

11()()(/2)4()

(),1,2,14m m m m m T h T h T h T h T h m -+-=??

?-=

=??- b.程序：

函数

function y=f(x) y=exp(x).*cos(x);

求解

M=1;h=3-1;err=1; R=zeros(4,4);j=0;eps=10^(-6); R(1,1)=h*(f(1)+f(3))/2; while (err>eps) j=j+1; h=h/2;

x=1+h:2*h:3-h;

R(j+1,1)=R(j,1)/2+h*sum(f(x)); for k=1:min(3,j)

R(j+1,k+1)=R(j+1,k)+(R(j+1,k)-R(j,k))/(4^k-1); end

if (j>3)

err=abs(R(j+1,4)-R(j,4)); end end

quad=R(j,4); disp(quad);

c.运行结果：

-10.4031

1、分别用Euler 方法（h=0.025）、改进的Euler 方法（h=0.05）、经典的R-K 方法(h=0.1）求初值问题的数值解计算0.1,0.2,0.3,0.4,0.5 处解的近似值。并与精确解比较，分析误差。

()'2100.5

y y x y =-+=

解：

(1)Euler 方法： a.原理：

欧拉公式：1

00(,)0,1,2,,1

()k k k k y y hf x y k n y x y +=+=-??=?

b.程序：

clc

x0=0;y0=0.5;h=0.025;n=5; x=x0;y=y0; for i=1:n

y=y+h*(y-x^2+1); x=x+4*h; disp('x=');disp(x); disp('y=');disp(y); end

c.运行结果：

(2) 改进的Euler 方法 a.原理：

改进的Euler 方法公式：()()

1211

12,,/2()n n n n n n k f x y k f x h y hk y y h k k +=??

=++??=++?

b.程序：

clc

x0=0;y0=0.5;h=0.05;n=5;

x=x0;y=y0;

for i=1:n

y1=y+h*(y-x^2+1);

y2=y+h*(y1-(x+2*h)^2+1);

y=(y1+y2)/2;

x=x+2*h;

disp('x=');disp(x);

disp('y=');disp(y);

end

c.运行结果：

(3)经典的R-K方法

a.原理：

经典的R-k公式：

1

21

32

43

11234

(,)

/2,

2

/2,

2

(,)

(22)

6

n n

n n

n n

n n

n n

k f x y

h

k f x h y k

h

k f x h y k

k f x h y hk

h

y y k k k k +

=

?

?

??

?=++

?

???

?

???

=++

? ?

??

?

?=++

?

?=++++

??

b.程序：

clc

x0=0;y0=0.5;h=0.1;n=5;

x=x0;y=y0;

for i=1:n

K1=h*(y-x^2+1);

K2=h*(y+K1/2-(x+h/2)^2+1);

K3=h*(y+K2/2-(x+h/2)^2+1);

K4=h*(y+K3-(x+h)^2+1);