数值分析课程设计报告.docx

数值分析课程设计报告.docx

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

数值分析课程设计报告

数值分析课程设计报告

院系理学院

班级11级应用数学

学生姓名汪佩祺

学号201130512276

任课教师雷秀仁

完成日期2013年8月23日

目录

设计题一

（1）问题……………………………………………………………………………………1

（2）设计思路………………………………………………………………………………1

（3）程序清单………………………………………………………………………………1

（4）程序运行与结果………………………………………………………………………5

（5）结果分析………………………………………………………………………………6

设计题二

（1）问题……………………………………………………………………………………7

（2）设计思路………………………………………………………………………………7

（3）程序清单………………………………………………………………………………7

（4）程序运行与结果………………………………………………………………………9

（5）结果分析………………………………………………………………………………11

设计题四

（1）问题……………………………………………………………………………………11

（2）设计思路………………………………………………………………………………12

（3）程序清单………………………………………………………………………………12

（4）程序运行与结果………………………………………………………………………12

（5）结果分析………………………………………………………………………………13

设计题五

（1）问题……………………………………………………………………………………14

（2）设计思路………………………………………………………………………………14

（3）程序清单………………………………………………………………………………14

（4）程序运行与结果………………………………………………………………………20

（5）结果分析………………………………………………………………………………24

心得体会…………………………………………………………………………………25

自我评价…………………………………………………………………………………26

设计题一

1）问题

编写解线性代数方程组的列主元高斯消去法的函数，并调用该函数计算某个9阶以上的非奇异阵A的逆矩阵。

通过计算AA-1检查答案，并与使用inv（A）所得结果和运行时间进行比较。

2）设计思路

分为选主元、消去、回代和矩阵求逆四个过程：

对于第一列，先选择最大的元素作为主元，并将该行于第一行对换，再进行消去过程。

……

对于第k列，对比由k到n的元素最大值作为主元，并将该行与第k行对换，再进行消去过程。

最终，该行列式转化为一个上三角阵，紧接着进行回代过程。

3）程序清单

M文件：

%Gaosixiaoqu.m

function[x]=Gaosixiaoqu（a,b）

n=length（a）;

x=zeros（n,1）;

a=[ab];

form=1:

n-1

max=m;

fori=m+1:

n%选取主元

ifa（i,m）>a（max,m）

max=i;

end

end

temp=a（m,m:

n+1）;

a（m,m:

n+1）=a（max,m:

n+1）;

a（max,m:

n+1）=temp;

fori=m+1:

n%消去过程

a（i,m）=-a（i,m）/a（m,m）;

a（i,m+1:

n+1）=a（i,m+1:

n+1）+a（i,m）*a（m,m+1:

n+1）;

end

end

x（n,1）=a（n,n+1）/a（n,n）;%回代过程

fori=n-1:

-1:

1

sum=0;

forj=i+1:

n

sum=sum+x（j,1）*a（i,j）;

end

x（i,1）=（a（i,n+1）-sum）/a（i,i）;

end

命令文件：

>>a=[132914876;234970181;275098924;321044627;120938541;223455120;512345673;119234512;234127413]

>>b=[2;5;0;3;1;2;3;4;8]

>>Gaosixiaoqu（a,b）

ans=

-14.6399

-5.9055

1.2273

-5.9630

2.4398

7.2287

-6.3115

10.6150

5.0669

>>b1=[1;0;0;0;0;0;0;0;0]

>>x1=Gaosixiaoqu（a,b1）

x1=

2.1040

1.0097

-0.0492

0.8453

-0.5371

-0.9420

0.6956

-1.3804

-0.5398

>>b2=[0;1;0;0;0;0;0;0;0]

>>x2=Gaosixiaoqu（a,b2）

x2=

-1.7759

-0.7142

0.0719

-0.6657

0.4042

0.7226

-0.6298

1.2118

0.5046

......

>>x=[x1:

x2:

x3:

x4:

x5:

x6:

x7:

x8:

x9]

ans=

2.1040-1.77590.1792-0.6625-1.99404.12950.3231-0.3950-1.7044

1.0097-0.71420.1930-0.4288-0.95211.73660.0459-0.3053-0.5631

-0.04920.0719-0.0321-0.0187-0.0043-0.1076-0.01540.09990.1110

0.8453-0.66570.0458-0.2217-0.68551.61360.0161-0.1118-0.7253

-0.53710.4042-0.01240.25410.4710-0.8117-0.07700.11390.2073

-0.94200.7226-0.12870.28180.8782-1.6985-0.08600.12960.8641

0.6956-0.62980.1619-0.2518-0.55031.19180.1342-0.0366-0.7359

-1.38041.2118-0.14400.34921.2259-2.7525-0.03780.18381.2409

-0.53980.5046-0.14440.36190.4777-1.0513-0.16050.08200.5395

>>inv（a）

ans=

2.1040-1.77590.1792-0.6625-1.99404.12950.3231-0.3950-1.7044

1.0097-0.71420.1930-0.4288-0.95211.73660.0459-0.3053-0.5631

-0.04920.0719-0.0321-0.0187-0.0043-0.1076-0.01540.09990.1110

0.8453-0.66570.0458-0.2217-0.68551.61360.0161-0.1118-0.7253

-0.53710.4042-0.01240.25410.4710-0.8117-0.07700.11390.2073

-0.94200.7226-0.12870.28180.8782-1.6985-0.08600.12960.8641

0.6956-0.62980.1619-0.2518-0.55031.19180.1342-0.0366-0.7359

-1.38041.2118-0.14400.34921.2259-2.7525-0.03780.18381.2409

-0.53980.5046-0.14440.36190.4777-1.0513-0.16050.08200.5395

验证得x=inv（a），二者结果一致。

4）程序运行与结果

M文件

过程图1

过程图2

结果图

5）结果分析

本程序利用列主元消去法解矩阵的逆主要是通过解AX=b，根据矩阵求逆

过程中的初等变换为依据，将单位矩阵分解成n个列向量，依次求出AXi=bi的解，从而得到inv（a）=[x1:

x2:

x3:

x4:

x5:

x6:

x7:

x8:

x9]，经过验证与结果完全一致。

cond（a）=62.4754，是个大数，故a，b的扰动引起的解x的相对误差就比较大，所以稳定性比较弱。

优点

这个程序过程清晰易懂，明确的分清了换行、消去过程、回代过程等，步骤明确，将高斯法解线性方程组运用到矩阵的求逆上，是对书本知识的一个应用之一。

缺点

但是，求逆的过程中没有使用较好的方法，浪费了时间。

设计题二

1）问题

对于迭代法

，它显然有不动点

。

试不用判定收敛阶的定理，设计1至2个数值实验（其中必须有一个不是直接用收敛阶的定义）得到收敛阶数的大概数值。

2）设计思路

I、直接用定义求收敛阶：

迭代过程

显然有不动点

。

记迭代误差

，如果存在实数p>=1和非零常数c，使得

，则迭代过程是p阶收敛的。

II、间接用定义求收敛阶：

通过分析题中迭代法，可以知道只需判断e[k+1]的值即可，因为e[k]在非收敛阶情况下，去极限值恒为0。

3）程序清单

m文件

I、直接用定义求收敛阶：

%Qiujie.m

functionn=Qiujie（x0）

m=0;

symsx;

x1=0.99*x-x^2;

k=（abs（x0-x1））/（abs（x0-x））^m;

whilelimit（k,x,x0）==0

m=m+1;

k=（abs（x0-x1））/（abs（x0-x））^m;

x=x1;

end

n=m;

II、间接用定义求收敛阶：

%Aqiujie.m

function[]=Aqiujie（）

x=0;

e=x^2-0.99*x;

ek=-x;

m=0;

whileabs（e）==0%判断条件

m=m+1;

e=diff（e）;

end

disp（'收敛阶数为p='）;%输出收敛阶和收敛的性质

disp（p）;

ifm==1

disp（'迭代法线性收敛'）;

elseifm~=1

disp（'迭代法超线性收敛'）

end

end

命令文件

方法一：

>>Qiujie（0）

ans=

1

方法二：

>>Aqiujie

收敛阶数为p=

1

迭代法线性收敛

4）程序运行与结果

方法一m文件

方法二m文件

方法一运行结果

方法二运行结果

5）结果分析

方法一：

该收敛公式对于x=0是1阶收敛的。

方法二：

从输出结果可以知道，这个迭代法在不动点x=0处是1阶收敛的。

优点

一、提供了两种满足题目要求的方法。

二、程序输出都有相对详细的提示语，可以比较清楚地看到输出的是几阶收敛。

缺点

方案具有明显的局限性，仅仅是围绕定义来设计。

设计题四

1）问题

某飞机头部的光滑外形曲线的型值点坐标由下表给出:

0

1

2

3

4

5

6

7

8

9

10

0

70

130

210

337

578

776

1012

1142

1462

1841

0

57

78

103

135

182

214

244

256

272

275

试建立其合适的模拟曲线（未必是用拟合方法），并求在点x＝100，250，400，500，800处的函数值y及一阶、二阶导数值y’，y”。

绘出模拟曲线的图形。

2）设计思路

先根据给定的数据用画图指令画出大致图像，分析其大概形状，然后建立一条比较合适的曲线去逼近原函数。

3）程序清单

formatlong

x1=[0701302103375787761012114214621841];

y1=[05778103135182214244256272275];

x0=[100250400500800];

plot（x1,y1）

polyfit（x1,y1,2）;%拟合出给出数据函数

y2=-0.0001*x1.^2+0.3245*x1+28.0193;

nihehanshu（x0）

plot（x1,y1,'b',x1,y2,'r'）

title（['模拟曲线和数据曲线']）

xlabel（'x'）;ylabel（'y'）;

legend（{'数据曲线','拟合曲线'}）;

gridon;

4）程序运行与结果

点x＝100，250，400，500，800处的函数值y及一阶、二阶导数值y’，y”如下：

x

100

250

400

500

800

y

59.4693

102.8943

141.8193

165.2693

223.6193

y’

0.1138

-0.2019

-0.5178

-0.7284

-1.3602

y”

-0.002

-0.0021

-0.0021

-0.0021

-0.002

5）结果分析

通过Matlab可算得一元二次函数的的表达式：

y=-0.0001*x^2+0.3245*x+28.0193。

优点

I、程序简单易懂。

II、经过对比，拟合结果曲线有比较好的拟合，结果有一定程度的可靠性。

缺点

程序过于简单，没有深入思考。

结果不一定是最佳的拟合曲线。

设计题五

1）问题

给定初值问题

其精确解为

，分别按下列方案求它在节点

处的数值解及误差。

比较各方法的优缺点，并将计算结果与精确解做比较（列表、画图）。

（方案I）欧拉法，步长h=0.025,h=0.1；

（方案II）改进的欧拉法，步长h=0.05,h=0.1；

（方案III）四阶经典龙格—库塔法，步长h=0.1。

2）设计思路

1、设计出题目要求的三种方案的m文件以及微分方程的m文件。

2、代入初值得到各方法的数值解。

3、对比精确解，进行误差分析。

4、作图进行对比。

3）程序清单

微分方程m文件

%fun.m

functionz=fun（x,y）

z=2/x*y+x.^2*exp（x）;

欧拉法m文件

%Euler.m

functionx=Euler（f,x0,y0,xn,M）

%f为一阶常微分方程的一般表达式的右端函数;x0，y0为初始条件;xn为取值范围的一个端点;M为区间的个数

x=zeros（1,M+1）;%x为Xn构成的向量

y=zeros（1,M+1）;%y为Yn构成的向量

x

（1）=x0;

y

（1）=y0;

h=（xn-x0）/M;

forn=1:

M

x（n+1）=x（n）+h;

y（n+1）=y（n）+h*feval（f,x（n）,y（n））;

end

T=[x',y']

改进欧拉法m文件

%GaijinEuler.m

functionx=GaijinEuler（f,x0,y0,xn,M）

x=zeros（1,M+1）;%x为Xn构成的向量

y=zeros（1,M+1）;%y为Yn构成的向量

x

（1）=x0;

y

（1）=y0;

h=（xn-x0）/M;

forn=1:

M;

x（n+1）=x（n）+h;

z0=y（n）+h*feval（f,x（n）,y（n））;

y（n+1）=y（n）+h/2*（feval（f,x（n）,y（n））+feval（f,x（n+1）,z0）;

end

T=[x',y']

四阶经典龙格—库塔法m文件

%FRonK.m

functionR=FRonK（f,x0,y0,xn,M）;

h=（y0-x0）/M;

x=zeros（1,M+1）;

y=zeros（1,M+1）;

T=x0:

h:

y0;

y

（1）=xn;

forn=1:

M

k1=h*feval（f,T（n）,y（n））;

k2=h*feval（f,T（n）+h/2,y（n）+k1/2）;

k3=h*feval（f,T（n）+h/2,y（n）+k2/2）;

k4=h*feval（f,T（n）+h,y（n）+k3）;

y（n+1）=y（n）+（k1+2*k2+2*k3+k4）/6;

end

R=[T',y']

命令文件

欧拉法

h=0.025时

>>clearall;

>>x0=1;y0=0;xn=2;M=40;

x=Euler（'fun',x0,y0,xn,M）

T=

1.00000

1.02500.0680

1.05000.1445

1.07500.2301

1.10000.3255

1.12500.4311

1.15000.5478

1.17500.6760

1.20000.8165

1.22500.9701

1.25001.1374

1.27501.3192

1.30001.5164

1.32501.7297

1.35001.9601

1.37502.2085

1.40002.4757

1.42502.7629

1.45003.0709

1.47503.4009

1.50003.7539

1.52504.1311

1.55004.5337

1.57504.9630

1.60005.4201

1.62505.9065

1.65006.4235

1.67506.9725

1.70007.5551

1.72508.1728

1.75008.8272

1.77509.5200

1.800010.2529

1.825011.0277

1.850011.8464

1.875012.7107

1.900013.6228

1.925014.5847

1.950015.5985

1.975016.6667

2.000017.7914

x=

Columns1through14

1.00001.02501.05001.07501.10001.12501.15001.17501.20001.22501.25001.27501.30001.3250

Columns15through28

1.35001.37501.40001.42501.45001.47501.50001.52501.55001.57501.60001.62501.65001.6750

Columns29through41

1.70001.72501.75001.77501.80001.82501.85001.87501.90001.92501.95001.97502.0000

h=0.1时

>>clearall;

x0=1;y0=0;xn=2;M=10;

x=Euler（'fun',x0,y0,xn,M）

T=

1.00000

1.10000.2718

1.20000.6848

1.30001.2770

1.40002.0935

1.50003.1874

1.60004.6208

1.70006.4664

1.80008.8091

1.900011.7480

2.000015.3982

x=

1.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.0000

改进欧拉法

h=0.05时

>>clearall;

x0=1;y0=0;xn=2;M=20;

x=GaijinEuler（'fun',x0,y0,xn,M）

T=

1.00000

1.05000.1532

1.10000.3449

1.15000.5801

1.20000.8643

1.25001.2032

1.30001.6031

1.35002.0710

1.40002.6142

1.45003.2406

1.50003.9589

1.55004.7785

1.60005.7092

1.65006.7621

1.70007.9487

1.75009.2816

1.800010.7744

1.850012.4418

1.900014.2993

1.950016.3641

2.000018.6542

h=0.1时

>>clearall;

x0=1;y0=0;xn=2;M=10;

x=GaijinEuler（'fun',x0,y0,xn,M）

T=

1.00000

1.10000.3424

1.20000.8583

1.30001.5927

1.40002.5983

1.50003.9364

1.60005.6789

1.70007.9092

1.800010.7245

1.900014.2374

2.000018.5789

四阶经典龙格—库塔法

h=0.1时

>>clearall;

x0=1;y0=0;xn=2;M=10;

R=FRonK（'fun',x0,y0,xn,M）

R=

1.00002.0000

0.90001.4105

0.80000.9647

0.70000.6348

0.60000.3975

0.50000.2327

0.40000.1239

0.30000.0570

0.20000.0205

0.10000.0049

0Inf

画图程序

精确解函数m文件

%f.m

functionz=f（x）;

z=x.^2*（exp（