数值分析课程设计报告Word文件下载.docx
《数值分析课程设计报告Word文件下载.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计报告Word文件下载.docx（28页珍藏版）》请在冰豆网上搜索。
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
temp=a（m,m:
n+1）;
a（m,m:
n+1）=a（max,m:
a（max,m:
n+1）=temp;
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:
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;
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]
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;
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]
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）
验证得x=inv（a），二者结果一致。
4）程序运行与结果
M文件
过程图1
过程图2
结果图
5）结果分析
本程序利用列主元消去法解矩阵的逆主要是通过解AX=b，根据矩阵求逆
过程中的初等变换为依据，将单位矩阵分解成n个列向量，依次求出AXi=bi的解，从而得到inv（a）=[x1:
x9]，经过验证与结果完全一致。
cond（a）=62.4754，是个大数，故a，b的扰动引起的解x的相对误差就比较大，所以稳定性比较弱。
优点
这个程序过程清晰易懂，明确的分清了换行、消去过程、回代过程等，步骤明确，将高斯法解线性方程组运用到矩阵的求逆上，是对书本知识的一个应用之一。
缺点
但是，求逆的过程中没有使用较好的方法，浪费了时间。
对于迭代法
，它显然有不动点
。
试不用判定收敛阶的定理，设计1至2个数值实验（其中必须有一个不是直接用收敛阶的定义）得到收敛阶数的大概数值。
2）设计思路
I、直接用定义求收敛阶：
迭代过程
显然有不动点
记迭代误差
，如果存在实数p>
=1和非零常数c，使得
，则迭代过程是p阶收敛的。
II、间接用定义求收敛阶：
通过分析题中迭代法，可以知道只需判断e[k+1]的值即可，因为e[k]在非收敛阶情况下，去极限值恒为0。
m文件
%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;
%Aqiujie.m
function[]=Aqiujie（）
x=0;
e=x^2-0.99*x;
ek=-x;
whileabs（e）==0%判断条件
e=diff（e）;
disp（'
收敛阶数为p='
）;
%输出收敛阶和收敛的性质
disp（p）;
ifm==1
disp（'
迭代法线性收敛'
elseifm~=1
迭代法超线性收敛'
）
命令文件
方法一：
Qiujie（0）
1
方法二：
Aqiujie
收敛阶数为p=
迭代法线性收敛
方法一m文件
方法二m文件
方法一运行结果
方法二运行结果
该收敛公式对于x=0是1阶收敛的。
从输出结果可以知道，这个迭代法在不动点x=0处是1阶收敛的。
一、提供了两种满足题目要求的方法。
二、程序输出都有相对详细的提示语，可以比较清楚地看到输出的是几阶收敛。
方案具有明显的局限性，仅仅是围绕定义来设计。
设计题四
某飞机头部的光滑外形曲线的型值点坐标由下表给出:
2
3
4
5
6
7
8
9
10
70
130
210
337
578
776
1012
1142
1462
1841
57
78
103
135
182
214
244
256
272
275
试建立其合适的模拟曲线（未必是用拟合方法），并求在点x＝100，250，400，500，800处的函数值y及一阶、二阶导数值y’，y”。
绘出模拟曲线的图形。
先根据给定的数据用画图指令画出大致图像，分析其大概形状，然后建立一条比较合适的曲线去逼近原函数。
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;
点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
通过Matlab可算得一元二次函数的的表达式：
y=-0.0001*x^2+0.3245*x+28.0193。
I、程序简单易懂。
II、经过对比，拟合结果曲线有比较好的拟合，结果有一定程度的可靠性。
程序过于简单，没有深入思考。
结果不一定是最佳的拟合曲线。
给定初值问题
其精确解为
，分别按下列方案求它在节点
处的数值解及误差。
比较各方法的优缺点，并将计算结果与精确解做比较（列表、画图）。
（方案I）欧拉法，步长h=0.025,h=0.1；
（方案II）改进的欧拉法，步长h=0.05,h=0.1；
（方案III）四阶经典龙格—库塔法，步长h=0.1。
1、设计出题目要求的三种方案的m文件以及微分方程的m文件。
2、代入初值得到各方法的数值解。
3、对比精确解，进行误差分析。
4、作图进行对比。
微分方程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））;
T=[x'
y'
]
改进欧拉法m文件
%GaijinEuler.m
functionx=GaijinEuler（f,x0,y0,xn,M）
%x为Xn构成的向量
%y为Yn构成的向量
M;
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）;
四阶经典龙格—库塔法m文件
%FRonK.m
functionR=FRonK（f,x0,y0,xn,M）;
h=（y0-x0）/M;
T=x0:
h:
y0;
y
（1）=xn;
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;
R=[T'
欧拉法
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时
M=10;
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
1.00001.10001.20001.30001.40001.50001.60001.70001.80001.90002.0000
改进欧拉法
h=0.05时
M=20;
x=GaijinEuler（'
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时
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
四阶经典龙格—库塔法
R=FRonK（'
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（
升级会员 