北京科技大学计算方法大作业

计算方法大作业

机械电子工程系 老师：廖福成

注：本文本只有程序题，证明题全部在手写已交到理化楼204了。

32. 证明方程 x?x?1?0

*|xk?1?xk|?10?3[1,2]x在上有一实根，并用二分法求这个根。要求。请给出程

序和运行结果。

1

证明： 设f(x)=x3-x-1

则f(1)= -1， f(2)= 5，f(1)*f(2)= -5<0 因此，方程在[1,2]上必有一实根。 二分法求解程序：

%预先定义homework2.m文件如下： function lc=homework2(x) lc=x^3-x-1;

在MALAB窗口运行： clear

a=1;b=2;tol=10^(-3);N=10000;

k=0; fa=homework2(a); % f 需事先定义 for k=1:N p=(a+b)/2; fp=homework2(p); if( fp==0 || (b-a)/2

if fa*fp<0 b=p; else a=p; end end k,p

程序运行结果： k = 10

p = 1.325195312500000

2

323. 用Newton迭代法求方程 x?2x?10x?20?0

的一个正根，计算结果精确到7位有效数字. 要求给出程序和运行结果. 解：

32'2x?1f(x)?x?2x?10x?20f(x)?3x?4x?10． 0取迭代初值 ,并设,则

f(x)x3?2x2?10x?20?(x)?x?'?x?2f(x)3x?4x?10 牛顿迭代函数为

xk3?2xk2?10xk?20xk?1?xk?23xk?4xk?10牛顿迭格式为：

Matlab程序如下： %定义zuoye3.m文件

function x=zuoye3(fname,dfname,x0,e,N) if nargin<5,N=500;end if nargin<4,e=1e-7;end x=x0;x0=x+2*e;k=0; while abs(x0-x)>e&k

x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); disp(x) end

if k==N,warning('已达上限次数');end 在Matlab窗口中执行：

3

zuoye3(inline('x^3+2*x^2+10*x-20'),inline('3*x^2+4*x+10'),1,1e-7) 结果如下： 1.41176470588235 1.36933647058824 1.36880818861753 1.36880810782137 ans =

1.36880810782137

3x?14. 用牛顿迭代法求方程x?x?1?0在0附近的根. 要求给出程序和运行结

果.

3'2f(x)?x?x?1f(x)?3x?1． 解：令：,则

f(x)x3?x?1?(x)?x?'?x?f(x)3x2?1 牛顿迭代函数为

xk3?xk?1xk?1?xk?23x?1 k牛顿迭格式为：

Matalb程序如下： %定义zuoye4.m文件

function x=zuoye4(fname,dfname,x0,e,N) if nargin<5,N=500;end if nargin<4,e=1e-7;end x=x0;x0=x+2*e;k=0;

4

while abs(x0-x)>e&k

x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); disp(x) end

if k==N,warning('已达上限次数');end

在Matlab窗口执行：

zuoye4(inline('x^3-x-1'),inline('3*x^2-1'),1,1e-7) 结果如下： 1.50000000000000 1.34782608695652 1.32520039895091 1.32471817399905 1.32471795724479 1.32471795724475 ans =

1.32471795724475

6. 编写用全主元Gauss消去法解线性方程组的程序，并求解

5

?0.02x1?x2?4x3?3x4?x5?11??x?x?2x?x?3x?1412345???4x1?2x2?3x3?3x4?x5?4??3x?x?3x?2x?4x?1612345???x1?3x2?x3?4x4?4x5?18

解：

Matlab程序如下：

A=[2 -1 4 -3 1;-1 1 2 1 3;4 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4] b=[11 14 4 16 18] function x=zuoye6(A,b) [n,v]=size(b);

D=[A,b;eye(n),zeros(n,v)],[s1,m]=size(D); for k=1:(n-1)

s=abs(A(k,k));p=k;q=k; for i=k:n for j=k:n if abs(A(i,j))>s

s=abs(A(i,j));p=i;q=j; end end end

if p>k t=D(k,:); D(k,:)=D(p,:); D(p,:)=t; end if q>k t1=D(:,k); D(:,k)=D(:,q); D(:,q)=t1; end h=D(k+1:n,k)/D(k,k);

D(k+1:n,k+1:m)=D(k+1:n,k+1:m)-h*D(k,k+1:m); D(k+1:n,k)=zeros(n-k,1);

6

end for k=n:-1:1

D(k,k:m)=D(k,k:m)/D(k,k);

for r=1:k-1 D(r,:)=D(r,:)-D(r,k)*D(k,:); end end

P=D(n+1:2*n,1:n);Q=D(1:n,n+1:m);x=P*Q 在Matlab窗口中执行：

A=[0.02 -1 4 -3 1;-1 1 2 1 3;4 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4]; b=[11 14 4 16 18]'; zuoye6(A,b) 运行结果如下： x =

2.94117647058824 -3.82352941176472 1.00000000000000 0.94117647058824 5.94117647058824

0?2?10??12?10??0?12?1?0?12??0?000?17. 用追赶法解线性方程组 ?

7

0??x1???x?0??2???x3??0????x4??1?2???x5??????????1??0??0?0??0 ?要求给出程序和运行结果. 解：

?1?1??00??2?100100000??2??0??0?130???100?00???2A???12?100???2??2?0?12?10???03100?????00?00?12?1?LU???????000?12?0?310????0?4??0???0??000?451???????00{Ly?b于是有求解Ax=b即为求解Ux?y，式中b=(1 0 0 0 0)T

??10000??1???1???1000????1???2??0?23100??y1????1???2???y2??0??1????y??3???0??00?310????y??40??3?1???4?????y5??????0????4??000?4??据

?51??1?? y=??5??

?5??2?1000?????1???6??03?1?????2?100????2???4??2?3?1???003?10???x???1??1?2?????x?2?0005?1????3?1???1??4??6??x3????x???4???3?4据??0000?5?????1???x5??=??5?? x=?1??6??

Matlab程序如下：

8

?43?10???05?4?1?006??5?? %定义zuoye7.m文件 function x=zuoye7(a,b,c,d) a1=[0;a];

n=length(b);q=zeros(n,1);p=zeros(n,1); %LU分解

q(1)=b(1);for k=2:n,p(k)=a1(k)/q(k-1); q(k)=b(k)-p(k)*c(k-1); end %解Ly=d

y=zeros(n,1);y(1)=d(1);for k=2:n, y(k)=d(k)-p(k)*y(k-1);end %解Ux=y

x=zeros(n,1); x(n)=y(n)/q(n);for k=n-1:-1:1,x(k)=(y(k)-c(k)*x(k+1))/q(k);end x

在Matlab窗口中执行： a=[-1 -1 -1 -1]'; b=[2 2 2 2 2]'; c=[-1 -1 -1 -1]'; d=[1 0 0 0 0]'; x=zuoye7(a,b,c,d) 运行结果如下： x = 0.83333333333333 0.66666666666667 0.50000000000000 0.33333333333333 0.16666666666667

9

9. 分别用Jacobi迭代法和Gauss-Seidel迭代法求解程组（编写程序）

?10x1+3x2?x3?14????2x1-10x2?3x3?-5? ?x?3x?10x?14?23?1?T取初值X（0）=（0,0,0），精确到小数后面四位。

解：

（1） Jacobi迭代法的Matlab程序：

% x0为初始向量，ep为精度，N为最大次数，x是近似解向量 A=[10 3 1;2 -10 3;1 3 10];

b=[14 -5 14];n=length(b);N=500;ep=1e-6;x0=zero(n,1); n=length(b);x0=zeros(n,1);x=zeros(n,1);k=0 while K

x(i)=(b(i)-A(I,[1:i-1,i+1:n])*x0(1:i-1,i+1:n))/A(i,i); end

if norm(x-x0,inf)

if k==N,Warning(‘已达到迭代次数上限’);end disp([‘k=’,num2str(k)]),x 运行结果如下： k=15 x=

1.000000327906423

10

0.999999801006905 1.000000327906432

（2）Gauss-Seidel迭代法的Matlab程序：

% x0为初始向量，ep为精度，N为最大次数，x是近似解向量 Format long;clear; A=[10 3 1;2 -10 3;1 3 10];

b=[14 -5 14];n=length(b);N=500;ep=1e-6;x0=zero(n,1);P=inf; %以下是Guass-Seidal迭代法程序

D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);dD=det(D); if dD==0

disp(‘请注意：因为对角矩阵D奇异，所以此方程组无解。’) else

disp(‘请注意：因为对角矩阵D非奇异，所以此方程组无解。’) iD=inv(D-L);B_GS=ID*U;d_GS=iD*b; end k=0; while k

if k==N,warning(‘已达到迭代次数上限’)；end disp([‘k=’,num2str(k)]),x,%D,U,L,B_GS,d_GS,%jx=A\\b, 运行结果如下：

11

请注意：因为对角矩阵D非奇异，所以此方程组有解。 k=9

x=1.000000043651052 1.000000031224555 0.999999986267528

10．编写幂法程序求矩阵

10.5??1?A??110.25???2??0.50.25?

?x1(??10?6)1按模最大的特征值和对应的特征向量。 解：幂法程序如下：

% A 为n 阶方阵，u0 为初始向量，

%epsilon 为上限，max1 为循环次数。lambda 返回按模最大的特征值，u 返回

对应的特征向量。

clear;format long;clc;max1=1000; epsilon=1e-6; A=[1 1 0.5;1 1 0.25;0.5 0.25 2]; u0=[1;1;1]; u=u0;v0=u0; lambda=0;k=0;err=1;R=[]; while ((kepsilon))

v=A*u; [m,j]=max(abs(v)); m=v(j); u=v/m; err=abs(lambda-m); lambda=m; k=k+1;

Temp=[k;m;v;u];R=[R Temp];

12

end

k, lambda, u, disp(R), xlswrite('d:\\\\1.xls', R,'sheet1', 'b1') ;

%[D,Lmd]=eig(A) %D 之列为A 的特征向量，Lmd 之对角线上为A 的特征值

运行结果如下： k = 23

lambda = 2.536527202038736 u = 0.748222175139906 0.649662222855908 1.000000000000000

11. 编写用原点位移加速反幂法程序求矩阵

? 7 3 -2??A?? 3 4 -1????-2 -1 3 ??

?10(??10) q?1.9最接近于的特征值和相应的特征向量。

解：

反幂法程序如下：

% A 为n 阶方阵，v0,u0 为初始向量，

%epsilon 为上限，max1 为循环次数，q 为近似特征值。lambda 为提高精度的

特征值，v 给出对应的特征向量。

13

clear;format long;max1=100; epsilon=1e-10;

A=[7 3 -2;3 4 -1;-2 -1 3]; u0=[1;1;1]; u=u0;v=u0;q=1.9; lambda=0;k=0;err=1;mu=0.5;

n=length(u0); z=zeros(n,1);A=A-q*eye(n); %LU 分解

n=length(u0);UU=zeros(n,n);L=eye(n,n); UU(1,:)=A(1,:); L(2:n,1)=A(2:n,1)/UU(1,1); for s=2:n

UU(s,s:n)=A(s,s:n)-L(s,1:s-1)*UU(1:s-1,s:n);

L(s+1:n,s)=(A(s+1:n,s)-L(s+1:n,1:s-1)*UU(1:s-1,s))/UU(s,s); end

while ((kepsilon)) [m,j]=max(abs(v)); m=v(j); u=v/m; %解下三角方程组Lz=u

z(1)=u(1);for j=2:n, z(j)=u(j)-L(j,1:j-1)*z(1:j-1); end %解上三角方程组Uv=z v(n)=z(n)/UU(n,n);

for j=n-1:-1:1, v(j)=(z(j)-UU(j,j+1:n)*v(j+1:n))/UU(j,j); end err=abs(1/m-1/mu); k=k+1; mu=m; end

disp(['迭代次数k= ',num2str(k)])

disp('要求的矩阵A 的特征值lambda='), lambda= q+1/m disp('要求的矩阵A 的特征向量u'), u

14

运行结果如下： 迭代次数k= 17

要求的矩阵A 的特征值lambda= lambda =

2.000000000009091

要求的矩阵A 的特征向量u u =

0.999999999973922 -0.999999999956935 1.000000000000000

12. 已知插值点

(-2.00,17.00), (0.00,1.00), (1.00,2.00), (2.00,17.00),

求三次插值多项式，并计算f(0.6). 解：

其Matlab程序如下： function yy=zuoye10(x,y,xx) m=length(x);n=length(y);

if m~=n, error('向量x与y的长度必须一致');end

15

s=0; for i=1:n

t=ones(1,length(xx)); for j=1:n if j~=i

t=t.*(xx-x(j))/(x(i)-x(j)); end end s=s+t*y(i); end yy=s;

在Matlab窗口中执行： x=[-2.00 0.00 1.00 2.00]; y=[17.00 1.00 2.00 17.00]; xx=0.6; zuoye10(x,y,xx) 结果如下： ans =

0.25600000000000

13．编制Newton插值法的通用Matlab程序，并求f(0.5)的近似值. 已知的数值如下表所示

xi00.20.40.60.8

f(xi)0.19950.39650.58810.7721.09461

16

解：

Newton插值法的通用Matlab程序：

x=[0 0.2 0.4 0.6 0.8]';y=[0.1995 0.3965 0.5881 0.7721 0.9461]';xx=0.5; m=length(x);n=length(y);

if m~=n, error('向量x与y的长度必须一致');end for j=2:n，for i=j:n

newton_chazhi(i,j)=(newton_chazhi(i-1,j-1)-newton_chazhi(i,j-1))/(x(1+i-j)-x(i)); end end

yy= newton_chazhi(1,1);

t=1;for i=2:n，t=t*(xx-x(i-1));yy=yy+ newton_chazhi(i,i)*t; end

%计算误差近似值

err=newton_chazhi(n,n);for i=1:n,err=err*(xx-x(i));end,err yy

运行结果为：

err =-2.343750000006213e-06 yy=0.681187500000000

??y??f(x,y)??y|x?x0?y020．编写解常微分方程?的四阶龙格库塔法通用程序.

解：

四阶龙格库塔法通用程序：

function [x,y]=zuoye20(dyfun,xspan,y0,h)

17

x=xspan(1):h:xspan(2); y(1)=y0;

for n=1:length(x)-1

k1=feval(dyfun,x(n),y(n));

k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6; end x=x'; y=y';

21. 利用上题的程序求解初值问题:

?dy2?2??xy,?dx3??y(0)?1x?[0,1.2]

32y?1?x的数值解,并将结果与准确解进行比较. 取h?0.4.

解：

Matlab程序：

function [x,y]=zuoye21(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0;

for n=1:length(x)-1

18

k1=feval(dyfun,x(n),y(n));

k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6; end x=x' y=y'

在Matlab窗口中执行： dyfun=inline('2*x*y^(-2)/3'); xspan=[0,1.2]; y0=1; h=0.4;

nark18(dyfun,xspan,y0,h) 结果如下： x =

0 0.40000000000000 0.80000000000000 1.20000000000000 y =1.00000000000000 1.05075062243354 1.17933175912031 1.34631542500313

19

计算值与准确值比较如下：

22. 编写四阶龙格库塔法程序，并求解洛伦兹型系统（不许调用现成的龙格库塔法程序软件包）：

????x??y??yz?x???rx?qy?sxz?y?z????bz??xy

解：为了便于求解，这里首先为各个参数赋予初值，实际中根据需要修改程序中的参数即可。

取??0.25，??0.06，b?0.4，??0.5，r?120，q?1.3，s?1.5，???20.取初始点为(0.005,0.4596,?0.1146). 程序如下：

%预定义homework22.m文件 function f=homework22(t,y)

20

sgm=0.25;tao=0.06;xgm=0.5;r=120;q=1.3;s=1.5;b=0.4;u=-20;

f=[-sgm*y(1)+tao*y(2)+xgm*y(2)*y(3);r*y(1)-q*y(2)+s*y(1)*y(3);-b*y(3)+u*y(1)*y(2)];

主程序为：

y0=[0.0050 0.4596 -0.1146];h=0.005;

x=0:h:1000;y=zeros(length(y0),length(x));y(:,1)=y0'; %四阶龙格库塔程序 for n=1:(length(x)-1)

k1=homework22(x(n),y(:,n));

k2=homework22 (x(n)+h/2,y(:,n)+h/2*k1); k3=homework22 (x(n)+h/2,y(:,n)+h/2*k2); k4=homework22(x(n+1),y(:,n)+h*k3); y(:,n+1)=y(:,n)+h*(k1+2*k2+2*k3+k4)/6; end x=x';y=y';

%绘制图像

figure(1);plot(y(:,1),y(:,2));xlabel('x');ylabel('y'); figure(2);plot(y(:,2),y(:,3));xlabel('y');ylabel('z'); figure(3);plot(y(:,1),y(:,3));xlabel('x');ylabel('z'); figure(4);plot3(y(:,1),y(:,2),y(:,3));xlabel('x'); ylabel('y'); zlabel('z');

21

程序运行结果为绘制的其函数x-y-x的三维图像如下图所示：

其在各个坐标平面上的投影，即x-y、y-z、x-z的关系图像分别如下所示：

22

本文来源：https://www.bwwdw.com/article/bgtw.html