收藏
下载资源下载资源 加入VIP,免费下载

数值分析实验报告2.docx

上传人：b****5 文档编号：6997583 上传时间：2023-01-15 格式：DOCX 页数：16 大小：33.07KB
下载 相关 举报
数值分析实验报告2.docx_第1页
第1页 / 共16页
数值分析实验报告2.docx_第2页
第2页 / 共16页
数值分析实验报告2.docx_第3页
第3页 / 共16页
数值分析实验报告2.docx_第4页
第4页 / 共16页
数值分析实验报告2.docx_第5页
第5页 / 共16页
点击查看更多>>
下载资源
资源描述

数值分析实验报告2.docx

《数值分析实验报告2.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告2.docx（16页珍藏版）》请在冰豆网上搜索。

数值分析实验报告2.docx

数值分析实验报告2

 

数值分析上机实验报告

（2）

（注：

本实验报告中所有程序均为MATLAB语言程序）

 

班级：

姓名：

学号：

 

六章

2、用比例求根法求

在区间（1，

/2）内的一个根，直到近似根

满足精度|

（

）|<0.00001时终止运算。

程序：

clear

formatlong

x=1;

f=inline（'1-x*sin（x）'）;

a=1;

b=pi/2;

whileabs（a-b）>0.00001

iff（（a+b）/2）<0

b=（a+b）/2;

x=b;

end

iff（（a+b）/2）>0

a=（a+b）/2;

x=a;

end

iff（（a+b）/2）==0

x=（a+b）/2;

break

end

end

x

运行结果：

x=

1.114153168596455

4、比较一下两种求

+10

-2=0的根到三位有效数字所需的计算量。

（1）在区间（0,1）内用二分法；

（2）用迭代法

=（2-

）/10,取初值

=0。

（1）程序：

formatlong

a=0;

b=1;

R=1;

k=0;

f=inline（'exp（x）+10*x-2'）;

whileR>1e-5

c=（a+b）/2;

iff（a）*f（c）>0;

a=c;

else

b=c;

end

R=b-a;

k=k+1;

end

k

x=c

运行结果：

k=

17

x=

0.09052276611328

首先初步估计精确值在0.09左右，故取误差限选择为10^-5，，得到的计算次数为17次。

（2）程序：

formatlong

f=inline（'（2-exp（x））/10'）;

disp（'x='）;

x=feval（f,0）;

disp（x）;

Eps=1E-5;

i=1;

while1

x0=x;

i=i+1;

x=feval（f,x）;

disp（x）;

ifx>1E10;

break;

end;

ifabs（x-x0）

break;

end

end

i,x

输出结果：

x=

0.10000000000000

0.08948290819244

.0906********

0.09051261667437

0.09052646805264

0.09052495168284

i=

6

x=

0.09052495168284

若要保留三位有效数字，则误差限应选择为10^-5，得到的迭代次数为5次。

7、用下列方法求

=0在

=2附近的根，根的准确值

=1.87938524…,要求计算结果准确到四位有效数字。

（1）用newton法；

（2）用弦截法，取

=2，

=1.9；

（3）用抛物线法，取

=1，

=3，

=2。

（1）Newton法程序：

clear

m=input（'m='）;

x

（1）=2;

f=inline（'x^3-3*x-1'）;

f1=inline（'3*x^2-3'）;

forn=1:

m

x（n+1）=x（n）-f（x（n））/f1（x（n））;

end

x（n+1）

运行结果：

m=10

ans=

1.879385241571817

（2）弦截法程序：

clear

m=input（'m='）;

x

（1）=2;

x

（2）=1.9;

f=inline（'x^3-3*x-1'）;

forn=2:

m

x（n+1）=x（n）-（f（x（n））/（f（x（n））-f（x（n-1））））*（x（n）-x（n-1））;

end

x（n+1）

运行结果：

m=9

ans=

1.879385241571817

（3）抛物线法程序：

首先建立M文件：

functionw=W（a,b）

f=inline（'x^3-3*x-1'）;

w=（f（b）-f（a））/（b-a）;

再有程序：

clear

m=input（'m='）;

x

（1）=1;

x

（2）=3;

x（3）=2;

f=inline（'x^3-3*x-1'）;

forn=3:

m

q（n）=W（x（n）,x（n-1））;

p（n）=（W（x（n）,x（n-2））-W（x（n）,x（n-1）））/（x（n-2）-x（n-1））;

w（n）=q（n）+p（n）*（x（n）-x（n-1））;

s（n）=sqrt（w（n）^2-4*f（x（n））*p（n））;

ifw（n）+s（n）>0

x（n+1）=x（n）-2*f（x（n））/（w（n）+sqrt（w（n）^2-4*f（x（n））*p（n）））;

end

ifw（n）-s（n）>0

x（n+1）=x（n）-2*f（x（n））/（w（n）-sqrt（w（n）^2-4*f（x（n））*p（n）））;

end

end

x=x（n+1）

运行结果：

m=10

x=

1.879385241570596

13．应用Newton法于方程

导出求

的迭代公式，并求出

的值

程序：

clear

m=input（'m='）;

x

（1）=11;

f=inline（'1-115/x^2'）;

f1=inline（'2*115*x^（-3）'）;

forn=1:

m

x（n+1）=x（n）-f（x（n））/f1（x（n））;

end

sqrta=x（n+1）

运行结果：

m=10

sqrta=

10.723805294763608

七、八、九章

求解下面线性方程的解：

A是一个主对角线为4，主对角线上下各三条对角线为1的100阶的方阵，b是一个1到100的一列矩阵。

利用不同种方法就解线性方程组。

Lu分解程序：

clear

formatshort

n=input（'n='）;

a=ones（n）;

b=diag（diag（a,-1）,-1）;

c=diag（diag（a,-2）,-2）;

d=diag（diag（a,1）,1）;

e=diag（diag（a,2）,2）;

g=diag（diag（a,3）,3）;

h=diag（diag（a,-3）,-3）;

f=4*eye（n）;

A=b+c+d+e+f+g+h;

b=[1:

100]';

[L,U]=lu（A）;

x=U\（L\b）;

x'

运行结果：

n=100

ans=

Columns1through11

0.01700.18550.32260.42390.49450.58860.69640.80500.90330.99931.0978

Columns12through22

1.19961.30081.40061.49981.59971.69991.80011.90012.00002.09992.2000

Columns23through33

2.30002.40002.50002.60002.70002.80002.90003.00003.10003.20003.3000

Columns34through44

3.40003.50003.60003.70003.80003.90004.00004.10004.20004.30004.4000

Columns45through55

4.50004.60004.70004.80004.90005.00005.10005.20005.30005.40005.5000

Columns56through66

5.60005.70005.80005.90006.00006.10006.20006.30006.40006.50006.6000

Columns67through77

6.69996.80006.90017.00027.09997.19967.29987.40067.50097.59957.6976

Columns78through88

7.79907.90398.00548.09648.18638.29498.42488.52928.57808.61588.7831

Columns89through99

9.04009.18158.92308.75859.153710.391110.27678.85456.469110.266214.0251

Column100

17.3099

Cholesky分解：

clear

formatshort

n=input（'n='）;

a=ones（n）;

b=diag（diag（a,-1）,-1）;

c=diag（diag（a,-2）,-2）;

d=diag（diag（a,1）,1）;

e=diag（diag（a,2）,2）;

g=diag（diag（a,3）,3）;

h=diag（diag（a,-3）,-3）;

f=4*eye（n）;

A=b+c+d+e+f+g+h;

b=[1:

100]';

R=chol（A）;

x=R\（R'\b）;

x'

运行结果：

n=100

ans=

Columns1through11

0.01700.18550.32260.42390.49450.58860.69640.80500.90330.99931.0978

Columns12through22

1.19961.30081.40061.49981.59971.69991.80011.90012.00002.09992.2000

Columns23through33

2.30002.40002.50002.60002.70002.80002.90003.00003.10003.20003.3000

Columns34through44

3.40003.50003.60003.70003.80003.90004.00004.10004.20004.30004.4000

Columns45through55

4.50004.60004.70004.80004.90005.00005.10005.20005.30005.40005.5000

Columns56through66

5.60005.70005.80005.90006.00006.10006.20006.30006.40006.50006.6000

Columns67through77

6.69996.80006.90017.00027.09997.19967.29987.40067.50097.59957.6976

Columns78through88

7.79907.90398.00548.09648.18638.29498.42488.52928.57808.61588.7831

Columns89through99

9.04009.18158.92308.75859.153710.391110.27678.85456.469110.266214.0251

Column100

17.3099

Jacobi方法迭代

函数M文件：

function[y,n]=jacobi（A,b,x0,eps）

ifnargin==3

eps=1.0e-6;

elseifnargin<3

error

return

end

D=diag（diag（A））;

L=-tril（A,-1）;

U=-triu（A,1）;

B=D\（L+U）;

f=D\b;

y=B*x0+f;

n=1;

whilenorm（y-x0）>=eps

x0=y;

y=B*x0+f;

n=n+1;

end

程序：

clear

formatshort

n=input（'n='）;

a=ones（n）;

b=diag（diag（a,-1）,-1）;

c=diag（diag（a,-2）,-2）;

d=diag（diag（a,1）,1）;

e=diag（diag（a,2）,2）;

g=diag（diag（a,3）,3）;

h=diag（diag（a,-3）,-3）;

f=4*eye（n）;

A=b+c+d+e+f+g+h;

b=[1:

n]';

m=diag（a）;

[x,n]=jacobi（A,b,m,1.0e-3）;

x'

n

运行结果：

因该方程的系数矩阵A不满足相关条件，故Jacobi方法迭代对该方程不适用

Gauss-Serdel迭代方法

函数M文件：

function[y,n]=gauseidel（A,b,x0,eps）

ifnargin==3

eps=1.0e-6;

elseifnargin<3

error

return

end

D=diag（diag（A））;

L=-tril（A,-1）;

U=-triu（A,1）;

G=（D-L）\U;

f=（D-L）\b;

y=G*x0+f;

n=1;

whilenorm（y-x0）>=eps

x0=y;

y=G*x0+f;

n=n+1;

end

命令程序：

clear

formatlong

n=input（'n='）;

a=ones（n）;

b=diag（diag（a,-1）,-1）;

c=diag（diag（a,-2）,-2）;

d=diag（diag（a,1）,1）;

e=diag（diag（a,2）,2）;

g=diag（diag（a,3）,3）;

h=diag（diag（a,-3）,-3）;

f=4*eye（n）;

A=b+c+d+e+f+g+h;

b=[1:

n]';

m=diag（a）;

[x,n]=gauseidel（A,b,m,1.0e-10）;

x'

n

运行结果：

n=100

ans=

Columns1through10

0.01700.18550.32260.42390.49450.58860.69640.80500.90330.9993

Columns11through20

1.09781.19961.30081.40061.49981.59971.69991.80011.90012.0000

Columns21through30

2.09992.20002.30002.40002.50002.60002.70002.80002.90003.0000

Columns31through40

3.10003.20003.30003.40003.50003.60003.70003.80003.90004.0000

Columns41through50

4.10004.20004.30004.40004.50004.60004.70004.80004.90005.0000

Columns51through60

5.10005.20005.30005.40005.50005.60005.70005.80005.90006.0000

Columns61through70

6.10006.20006.30006.40006.50006.60006.69996.80006.90017.0002

Columns71through80

7.09997.19967.29987.40067.50097.59957.69767.79907.90398.0054

Columns81through90

8.09648.18638.29498.42488.52928.57808.61588.78319.04009.1815

Columns91through100

8.92308.75859.153710.391110.27678.85456.469110.266214.025117.3099

 

n=

46

 

幂法：

建立matlab的函数文件eig_power.m如下：

function[V,D]=eig_power（A）

Maxtime=100;

Eps=1E-5;

n=length（A）;

V=ones（n,1）;

k=0;

m0=0;

whilek<=Maxtime

v=A*V;

[vmax,i]=max（abs（v））;

m=v（i）;

V=v/m;

ifabs（m-m0）

break;

end

m0=m;

k=k+1;

end

D=m;

在命令窗口中调用函数文件eig_power.m：

a=linspace（4,4,100）;

c=linspace（1,1,99）;

d=linspace（1,1,98）;

e=linspace（1,1,97）;

f=diag（a）;

g=diag（c,1）;

h=diag（c,-1）;

i=diag（d,2）;

j=diag（d,-2）;

k=diag（e,3）;

l=diag（e,-3）;

A=f+g+h+i+j+k+l;

[V,D]=eig_power（A）;

V'

D

输出结果：

ans=

Columns1through9

0.55000.68000.81000.94000.97000.99001.00001.00001.0000

Columns10through18

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns19through27

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns28through36

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns37through45

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns46through54

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns55through63

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns64through72

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns73through81

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns82through90

1.00001.00001.00001.00001.00001.00001.00001.00001.0000

Columns91through99

1.00001.00001.00001.00000.99000.97000.94000.81000.6800

Column100

0.5500

D=

10

其中D是最大特征值，V则是特征向量。

 

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置：首页 > 高等教育 > 院校资料

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1