数值计算方法上机答案.docx

数值计算方法上机答案.docx

《数值计算方法上机答案.docx》由会员分享，可在线阅读，更多相关《数值计算方法上机答案.docx（34页珍藏版）》请在冰点文库上搜索。

数值计算方法上机答案

本科实验报告

课程名称：

实验项目：

实验地点：

机房

专业班级：

采矿1206学号：

**********

*******

指导教师：

2014年7月3日

一、求非线性方程的根。

1、求方程

在

附近的是根，要求精度满足

.（牛顿切线法）

f=inline（'x-cos（x）'）;%f（x）

df=inline（'1+sin（x）'）;%f'（x）

n=1;

x0=input（'x0='）;

del=input（'del='）;

N=input（'N='）;

fprintf（'\nkx（k）'）;

fprintf（'\n%2d%f',0,x0）;

F0=f（x0）;dF0=df（x0）;

whilen

ifdF0==0

fprintf（'导数为0,迭代无法继续进行.'）;

return;

end

x1=x0-F0/dF0;

F1=f（x1）;

dF1=df（x1）;

if（（abs（x1-x0）

fprintf（'\n\n结果：

%f\n',x1）;

return;

end

fprintf（'\n%2d%f',n,x1）;

n=n+1;

x0=x1;

F0=F1;

dF0=dF1;

end

fprintf（'\n\n%d次迭代后未达到精度要求.\n',N）;

NewtonIteration

x0=1.5

del=1e-4

N=100

kx（k）

01.500000

10.784472

20.739519

结果：

0.739085

2、求方程

在

附近的是根，求出具有思维有效数字的根近似值..（简单迭代法）

clear

clc

phi=inline（'（0.8+x^2）^（1/3）'）;%迭代函数

x0=input（'x0='）;

del=input（'del='）;

N=input（'N='）;n=1;

fprintf（'\n%2d%f',0,x0）;

whilen

x=phi（x0）;

ifabs（x-x0）

fprintf（'\n\n近似解=%f\n',x）;

return

end

fprintf（'\n%2d%f',n,x）;

n=n+1;

x0=x;

end

fprintf（'\n\n%fd次迭代后未达到精度要求.\n',N）

x0=1

del=1e-4

N=100

01.000000

11.216440

21.316116

31.363004

41.385180

51.395688

61.400671

71.403034

81.404155

91.404687

101.404939

111.405059

近似解=1.405116

100.000000d次迭代后未达到精度要求.

二、求解线性方程组（直接法或迭代法）

1、

（列主元素消元法）

a=input（'a='）%[2,2,1,-3,8;-2,1,-1,-3,1;8,-1,3,8,-1;10,4,4,3,8];

[p,n]=size（a）;

forw=1:

p

[x,y]=find（a（w:

p,w）==max（max（a（w:

p,w））））;

q=a（w,:

）;

a（w,:

）=a（x,:

）;

a（x,:

）=q;

end

forj=1:

（p-1）

fori=（j+1）:

p

a（i,:

）=a（j,j）/a（i,j）.*a（i,:

）-a（j,:

）;

end

end

m=p;

whilem>1

s（m）=a（m,n）;

j=p;

while（（j>2）&（j>=m+1）&（j

s（m）=s（m）-a（m,j）*x（j）;

j=j-1;

end

x（m）=s（m）/a（m,m）;

m=m-1;

x

end

a=[221-38;-21-1-31;8-138-1;104438]

a=

221-38

-21-1-31

8-138-1

104438

x=

1.0000-1.00002.0000-2.0000

function[x,det,flag]=Gauss（A,b）

[n,m]=size（A）;

nb=length（b）;

ifn~=m

error（'therowsandcloumsofAmustbeequal!

'）;

return;

end

ifm~=nb

error（'therowsandcloumsofAmustbeequalthelengthofb!

'）;

return;

end

flag='OK';det=1;x=zeros（n,1）;

fork=1:

（n-1）

max1=0;

fori=k:

n

ifabs（A（i,k））>max1

max1=abs（A（i,k））;r=i;

end

end

ifmax1<1e-6

flag='failure';return;

end

ifr>k

forj=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

fori=k+1:

n

m=A（i,k）/A（k,k）;

forj=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）;

ifabs（A（n,n））<1e-6

flag='faliure';return;

end

fork=n:

-1:

1

forj=k+1:

n

b（k）=b（k）-A（k,j）*x（j）;

end

x（k）=b（k）/A（k,k）;

end

x（k）=b（k）/A（k,k）;

end

X

Det

flag

A=vpa（[17.031-0.615-2.9911.007-1.0060;-134.211-1-2.16.3-1.7;00.513-0.51-1.5;4.5013.11-3.907-61.70512.178.999;0.101-0.812-0.017-0.914.9180.1;1234.5521.80],6）

b=vpa（[0.23-22.32254240.23629.304-117.818],6）;

A=

[17.031,-0.615,-2.991,1.007,-1.006,0]

[-1.0,34.211,-1.0,-2.1,6.3,-1.7]

[0,0.5,13.0,-0.5,1.0,-1.5]

[4.501,3.11,-3.907,-61.705,12.17,8.999]

[0.101,-0.812,-0.017,-0.91,4.918,0.1]

[1.0,2.0,3.0,4.5,5.0,21.8]

>>Gauss（A,b）

A=

[17.031,-0.615,-2.991,1.007,-1.006,0]

[-1.0,34.211,-1.0,-2.1,6.3,-1.7]

[0,0.5,13.0,-0.5,1.0,-1.5]

[4.501,3.11,-3.907,-61.705,12.17,8.999]

[0.101,-0.812,-0.017,-0.91,4.918,0.1]

[1.0,2.0,3.0,4.5,5.0,21.8]

ans=

1.0000

-2.0000

2.9999

-4.0001

5.0000

-6.0008

clear

clc

n=input（'n='）;

A=input（'A='）;

b=input（'b='）;

x=input（'x='）;

epsilon=input（'\n精度='）;

N=input（'\n最大迭代次数N='）;

fprintf（'\n%d:

',0）;

fori=1:

n

fprintf（'%f',x（i））;

end

%以下是迭代过程

fork=1:

N

%这是第k步迭代，迭代前的向量在x0[]中，迭代后的向量在x[]中；

normal=0;

fori=1:

n

t=x（i）;

x（i）=b（i）;

forj=1:

n

ifj~=i

x（i）=x（i）-A（i,j）*x（j）;

end

end

x（i）=x（i）/A（i,i）;

temp=abs（x（i）-t）;%求范数于迭代在同一个循环中；

iftemp>normal

normal=temp;%这里用的是无穷范数

end

end%第i不迭代结束；

fprintf（'\n%d:

',k）;

fori=1:

n

fprintf（'%f',x（i））;%输出迭代过程

end

ifnormal

return

end

end

fprintf（'\n\n迭代%d次后仍未求得满足精度的解\n',N）;

n=8

A=[17,1,4,3,-1,2,3,-7;210-17-211-4;-11-82-52-11;241-11134-1;1317-151-24;-217-1212-18;345128-192;5111-11-710]

b=[71;43;-11;-37;-61;52;-73;21]

x0=[0;0;0;0;0;0;0;0]

精度=1e-4

最大迭代次数N=100

0:

0.0000000.0000000.0000000.0000000.0000000.0000000.0000000.000000

1:

4.1764714.3000001.3750003.3636364.0666674.3333333.8421052.100000

2:

2.9225162.0834500.5552548.5693667.2030702.3916968.4193711.770708

3:

1.796703-1.159744-1.3226198.9108208.2232983.3576457.2909535.892632

4:

4.2118430.507676-1.2417916.9272028.8477911.5521346.9022046.149035

5:

4.8681201.868149-2.5922027.3796888.4218241.2958416.8839584.935878

6:

4.5347650.742110-2.3382647.8659948.5206412.9955406.6619884.531746

7:

4.1340460.204098-1.9753807.8475618.5522223.1606677.1482374.421535

8:

3.9348030.230528-2.0501637.8467958.3502823.0342277.1274764.968270

9:

4.1827670.456461-1.8566067.7032878.4770762.7098147.0649175.050721

10:

4.2398730.623716-2.0374607.7409308.4934282.5261337.0804484.900471

11:

4.2239090.527844-2.0911437.7663918.4818132.7290846.9876304.900388

12:

4.2294720.494477-2.0251787.7442958.4957602.7621067.0363614.834351

13:

4.1809570.483627-2.0503767.7731678.4616392.7712547.0548314.862722

14:

4.1877690.472362-2.0135997.7616418.4737412.7694887.0430224.896300

15:

4.1986420.499955-2.0210727.7553988.4794262.7237627.0548634.884616

16:

4.1986920.501255-2.0384687.7601448.4745152.7349437.0402094.891221

17:

4.2055730.498191-2.0286037.7557178.4802942.7405807.0410484.880464

18:

4.1993120.497105-2.0346057.7600978.4761272.7421377.0457024.877388

19:

4.1974990.492006-2.0308367.7600368.4757992.7481837.0430294.883474

20:

4.1991700.494820-2.0284517.7582898.4772642.7418767.0458104.882012

21:

4.1984870.496008-2.0317677.7593908.4761002.7413467.0445474.883555

22:

4.1998550.495605-2.0303137.7585468.4771302.7422207.0436914.883052

23:

4.1995870.496071-2.0311267.7588478.4768822.7416557.0445844.881763

24:

4.1990620.495235-2.0312187.7591598.4765432.7430467.0440424.882558

25:

4.1993170.495278-2.0304527.7588308.4768572.7425887.0444094.882330

26:

4.1991060.495514-2.0309467.7590208.4766402.7422837.0444594.882489

27:

4.1992540.495419-2.0307707.7589248.4767442.7424677.0442214.882645

28:

4.1993260.495564-2.0307947.7588918.4767882.7422487.0443704.882398

29:

4.1992300.495486-2.0309247.7589758.4767042.7424297.0442904.882484

30:

4.1992730.495441-2.0307827.7589238.4767582.7424537.0443054.882462

31:

4.1992410.495482-2.0308407.7589468.4767312.7423837.0443504.882449>>x

x=

4.1992

0.4955

-2.0308

7.7589

8.4767

2.7424

7.0444

4.8824

function[x,det,flag]=Gauss（A,b）

[n,m]=size（A）;

nb=length（b）;

ifn~=m

error（'therowsandcloumsofAmustbeequal!

'）;

return;

end

ifm~=nb

error（'therowsandcloumsofAmustbeequalthelengthofb!

'）;

return;

end

flag='OK';det=1;x=zeros（n,1）;

fork=1:

（n-1）

max1=0;

fori=k:

n

ifabs（A（i,k））>max1

max1=abs（A（i,k））;r=i;

end

end

ifmax1<1e-6

flag='failure';return;

end

ifr>k

forj=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

fori=k+1:

n

m=A（i,k）/A（k,k）;

forj=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）;

ifabs（A（n,n））<1e-6

flag='faliure';return;

end

fork=n:

-1:

1

forj=k+1:

n

b（k）=b（k）-A（k,j）*x（j）;

end

x（k）=b（k）/A（k,k）;

end

x（k）=b（k）/A（k,k）;

end

A=[11111;12345;1361015;14102035;15153570];b=[10000]

b=

10000

>>Gauss（A,b）

ans=

5

-10

10

-5

1

三、用追赶法求解下列方程组（6位有效数字）。

解1

clear

clc

n=input（'n='）;

a=input（'a='）;

b=input（'b='）;

c=input（'c='）;

%a

（1）=0;c

（1）=0;

x=input（'x='）;

s=zeros（n,1）;

t=s;

temp=0;

fork=1:

n

s（k）=b（k）-a（k）*temp;%temp

t（k）=c（k）/s（k）;

temp=t（k）;

end

temp=0;

fork=1:

n

x（k）=（x（k）-a（k）*temp）/s（k）;

temp=x（k）;

end

fork=n-1:

-1:

1

x（k）=x（k）-t（k）*x（k+1）;

end

fprintf（'\nx=\n'）;

fork=1:

n

fprintf（'%f\n',x（k））;

end

n=5

a=[0-1-1-1-1]

b=[44444]

c=[-1-1-1-10]

x=[100000200]

x=

27.051282

8.205128

5.769231

14.871795

53.717949

解2

clear

clc

n=input（'n='）;

a=input（'a='）;

b=input（'b='）;

c=input（'c='）;

%a

（1）=0;c

（1）=0;

x=input（'x='）;

s=zeros（n,1）;

t=s;

temp=0;

fork=1:

n

s（k）=b（k）-a（k）*temp;%temp

t（k）=c（k）/s（k）;

temp=t（k）;

end

temp=0;

fork=1:

n

x（k）=（x（k）-a（k）*temp）/s（k）;

temp=x（k）;

end

fork=n-1:

-1:

1

x（k）=x（k）-t（k）*x（k+1）;

end

fprintf（'\nx=\n'）;

fork=1:

n

fprintf（'%f\n',x（k））;

end

n=10

a=[0111111111]

b=[-2-2-2-2-2-2-2-2-2-2]

c=[1111111110]

x=[-0.5-1.5-1.5-1.5-1.5-1.5-1.5-1.5-1.5-0.5]

x=

6.500000

12.500000

17.000000

20.000000

21.500000

21.500000

20.000000

17.000000

12.500000

6.500000

四、求解下列积分（E=1E-6）

解1

%复化的梯形公式和复化色Simpson公式计算定积分

f=inline（'（sin（x）/x）'）;%函数表达式可以更换；

a=input（'a='）;%积分区间左端点

b=input（'b='）;%积分区间右端点

n=input（'n='）;%区间n等分；

h=（b-a）/n;

temp=f（a）;

xk=a;

fori=1:

n-1

xk=xk+h;

temp=temp+2*f（xk）;

end

temp=temp+f（b）;

temp=temp*h/2;

fprintf（'\n复化梯形公式计算结果：

%f',temp）;

fuhuatx

a=1

b=2

n=1000

复化梯形公式计算结果：

0.659330

解2

%复化的梯形公式和复化色Simpson公式计算定积分

f=inline（'exp（-x^2）'）;%函数表达式可以更换；

a=input（'a='）;%积分区间左端点

b=input（'b='）;%积分区间右端点

n=input（'n='）;%区间n等分；

h=（b-a）/n;

temp=f（a）;

xk=a;

fori=1:

n-1

xk=xk+h;

temp=temp+2*f（xk）;

end

temp=temp+f（b）;

temp=temp*h/2;

fprintf（'\n复化梯形公式计算结果：

%f',temp）;

fuhuatx

a=0

b=1

n=1000

复化梯形公式计算结果：

0.746824

五、求下列矩阵所有最大和最小特征值和相应特征向量。

解1

（幂法）

function[m,u,flag]=Pow（A,delta,max1）

ifnargin<3

max1=100

end

ifnargin<2

delta=1e-5;

end