数值分析课程课程设计Word格式文档下载.docx
《数值分析课程课程设计Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数值分析课程课程设计Word格式文档下载.docx（24页珍藏版）》请在冰豆网上搜索。
解该问题主要使用递推算法，关于椰子数目的变化规律可以设起初的椰子数为po,第一至五次猴子在夜里藏椰子后，椰子的数目分别为po,p1,p2,P3,P4再设最后每个人分得x个椰子，由题：
r17（pk1）（k=0,1,2,3,4）x-（p51）
55
所以P55x1,p<
p<
11利用逆向递推方法求解
R5妇1,（k=0,1,2,3,4）
4
输出结果
•
1023
15621
MATLAB码：
n=input（'
n='
）;
n=15621forx=1:
np=5*x+1;
fork=1:
5p=5*p/4+1;
end
结果分析：
此题的解题思想是在迭代法中，判断p为整数时，输出与p
ifp==fix（p）,breakendenddisp（[x,p]）
dn
1x
1.2设，Indx
05x
（1）从10尽可能精确的近似值出发，利用递推公式：
1In5In1—（n1,2,L20）n
计算机从Ii到I20的近似值；
（2）从I30较粗糙的估计值出发，用递推公式：
Ini-In—（n30,29,L,3,2）
55n
计算从Il到I20的近似值；
解：
首先第一步，估计Io和I30的值：
symsxn;
int（xA0/（5+x）,0,1）
ans=log
（2）+log（3）-log（5）
eval（ans）
ans=
0.1823
则取I。
为0.18
int（xA30/（5+x）,0,1）ans=
931322574615478515625*log
（2）+931322574615478515625*log（3）-931
322574615478515625*log（5）-79095966183067699902965545527073/46
5817912560
ans=
MATLAB中小数点后保留四位，由上面计算知道积分的值不为了零。
所以I30的取值为0.00001-0.0001
i=input（'
i='
i=0.18;
>
ifi>
=0.1&
&
i<
=0.2
forn=1:
1:
20
i=-5*i+1/n
elseifi>
0&
=0.0001
forn=30:
-1:
2
i=（-1/5）*i+1/（5*n）
i=
1.1336e+005
-5.6679e+005
2.8339e+006
-1.4170e+007
7.0848e+007
-3.5424e+008
1.7712e+009
-8.8560e+009
4.4280e+010
-2.2140e+011
同理输入积分初始值i=0时可以得i=0.0884
第二种方法所得的结果相对来说比较精确一些，也比较可靠因为第一种方法每一迭代都将最初的误差放大了五倍，使得最终的误差
越来越大；
而第一种方法经过每一次迭代都将误差缩小为初始误差的五分之一，使得最终的误差越来越小，因此相对来说比较可靠，性能较好。
1.3绘制Koch分形曲线
问题描述：
从一条直线段开始，将线段中间的三分之一部分用一个等边
三角形的另两条边代替，形成具有5个结点C0SJ的国写（图1-4）;
在新的图形中，又将图中每一直线段中间的询之丁部分都用一个等边三角形的另两条边代替，再次形成新的图形（图sit5）,cos时，图形中共有17
个结点。
这种迭代继续进行下去可以形成Koch分形曲线。
问题分析：
考虑由直线段（2个点）产生第一个图形（5个点）的过程，设p1和R分别为原始直线段的两个端点。
现在需要在直线段的中间依次插入三个点P2，R，P4产生第一次迭代的图形（图1-4）。
显然，只位于P点右端直线段的三分之一处，已点绕p2旋转60度（逆时针方向）而得到的，故可以处理为向量巳已经正交变换而得到向量P2P3,形成算法如下：
（1）
P2
P（任P）/3.
（2）
R
P2（P5R）/3.
（3）
P3
%（EE）AT.
在算法的第三步中，A为正交矩阵。
这一算法将根据初始数据（P和R点的坐标），产生图1-4中5个结点的
坐标。
这5个结点的坐标数组，组成一个5X2矩阵。
这一矩阵的第一行为为R的坐标，第二行为R的坐标，第二行为P2的坐标……第五行为P5的坐标。
矩阵的第一列元素分别为5个结点的x坐标，第二列元素分别为5个结点的y坐标。
问题思考与实验：
（1）考虑在Koch分形曲线的形成过程中结点数目的变化规律。
设第k次迭代产生结点数为nk,第k1迭代产生结点数为nk1,试写出nk和nk1之间的递推关系式；
（2）参考问题分析中的算法，考虑图1-4到图1-5的过程，即由第一次迭代的5个结点的结点坐标数组，产生第二次迭代的17个结点的结点坐标数组的算法；
（3）考虑由第k次迭代的nk个结点的结点坐标数组，产生第k1次迭代的nk1个结点的结点坐标数组的算法；
（4）设计算法用计算机绘制出如下的
（1）nk14nk3
Koch分形曲线（图1-6）
（2）（3辩法及（4）代码分析：
p=[00;
100];
%P为初始两个点的坐标，第一列为x坐标，第二列为y坐标
n=2;
%n为结点数
A=[cos（pi/3）-sin（pi/3）;
sin（pi/3）cos（pi/3）];
%旋转矩阵
fork=1:
d=diff（p）/3;
%diff%
m=4*n-3;
%
q=p（1:
n-1,:
p（5:
4:
m,:
）=p（2:
n,:
%p（2:
）=q+d;
p（3:
）=q+d+d*A'
;
%p（4:
）=q+2*d;
%n=m;
plot（p（:
1）,p（:
2））%
axis（[01003.5]）
计算相邻两个点的坐标之差，得到相邻两点确定的向量
则d就计算出每个向量长度的三分之一，与题中将线段三等分对应
迭代公式
以原点为起点，前n-1个点的坐标为终点形成向量
迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标
用向量方法计算迭代后处于
4k+2位置上的点的坐标
4k+3位置上的点的坐标
4k位置上的点的坐标
迭代后新的结点数目
2.1用高斯消元法的消元过程作矩阵分解。
设
2023
A181
2315
消元过程可将矩阵A化为上三角矩阵U,试求出消元过程所用的乘数
m21、m31、伉1并以如下格式构造下三角矩阵L和上三角矩阵U
12023
Lm211,Ua22）a23）
4
（2）
m31m321a33
验证：
矩阵A可以分解为L和U的乘积，即A=LU。
矩阵LU分解MATLAB代码：
functionhl=zhjLU（A）
[n,n]=size（A）;
RA=rank（A）;
ifRA~=n
disp（'
因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的
秩RA如下：
’）；
RA,hl=det（A）;
return
ifRA==n
forp=1：
n
h（p）=det（A（1：
p,1：
p））;
hl=h（1：
n）;
fori=1：
ifh（1,i）==0
因为A的各阶主子式等不等于零，所以A能进彳fLU
分解.A的秩RA和各阶顺序主子式值hl庆次如下：
RA,hl
ifh（1,i）~=0
因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl如下：
forj=1：
U（1,j）=A（1,j）;
fork=2:
fori=2:
forj=2:
L（1,1）=1;
L（i,i）=1;
j
L（2,1）=A（2,1）/U（1,1）;
L（i,1）=A（i,1）/U（1,1）;
L（i,k）=（A（i,k）-L（i,1:
k-1）*U（1:
k-1,k））/U（k,k）;
else
U（k,j）=A（k,j）-L（k,1:
k-1,j）;
RA,hl,U,L
仅上上代码保存为M文件，并在命令窗口输入
A=[2023;
181;
2-315];
b=[000]'
h=zhjLU（A）
程序运行结果:
L=
U=
1.0000
0020.00002.00003.0000
0.0500
1.0000007.90000.8500
0.1000
-0.40511.00000015.0443
实验验证：
可以直接使用MATLAB内置LU分解
[LU]=lu（A）输出结果与上程序输出结果一致。
2.2用矩阵分解方法求上题中A的逆矩阵。
记
1
b0,b2
1,b3
1分别求解方程组AXb1，AXb2，AX
由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。
对第一
个方程组，由于A=LU，所以先求解下三角方程组LY知再求解上三角方程组UXY,则可得逆矩阵的第一列列向量；
类似可解第二、第三方程组，得逆矩阵的第二列列向量的第三列列向量。
由三个列向量拼装可得逆矩阵A1。
解:
MATLAB代码如下:
b1=[1;
0;
0];
b2=[0;
1;
b3=[0;
1];
A=[20,2,3;
1,8,1;
2,-3,15];
L=[1,0,0;
0.05,1,0;
0.1,-0.4051,1];
U=[2023;
07.90.85;
0015.0443];
Y1=L\b1
X1=U\Y1
Y2=L\b2X2=U\Y2Y3=L\b3
X3=U\Y3
Y1=
-0.0500
-0.1203
X1=
0.0517
-0.0055
-0.0080
Y2=
[X1X2X3]ans=
0.0517-0.0164-0.0257
-0.00550.12370.1165
-0.00800.02690.0934
而：
inv（A）=[X1X2X3]得证
0.4051
X2=
-0.0164
0.1237
0.0269
Y3=
1.4051
X3=-0.02570.11650.0934
2.3验证希尔伯特矩阵的病态性：
对于三阶矩阵
11/21/3
H1/21/31/4
1/31/41/5
取右端向量b【11/613/1247/6°
]t,验证：
（1）向量X[X1X2X3]T[111]T是方程组HXb的准确解；
（2）取右端向量b的三位有效数字得b[1.83L080.783]T,求方程组的准确解X,并与X的数据I1,1,1】'
作比较。
说明矩阵的病态性。
⑴
H=[11/21/3;
1/21/31/4;
1/31/41/5];
X=[1；
1；
1]；
b=H*X
b=
1.8333
1.0833
0.7833
与题中相同
（2）先求出解X,与数据[〔，"
亍作比较
b=[1.83;
1.08;
0.783];
X=H\b
X=
1.0800
0.5400
1.4400
与[1,1,1]相差较大，矩阵为病态矩阵
3.
1用泰勒级数的有限项逼近正弦函数
sinx,x[0,]
x,x[0,/2]
3——
xx/6,x[0,/2]xx3/6x5/120,x
用计算机绘出上面四个函数的图形。
MATLAB代码如下
（1）symsx;
taylor（sin（x））x=0:
0.01*pi:
piplot（x,sin（x））
⑵
symsx;
taylor（x）
x=0:
pi/2
plot（x）
⑶
taylor（x-xA3/6）
fplot（'
x-xA3/6'
[0pi/2]）（4）
taylor（x-xA3/6+xA5/120）
x-xA3/6+xA5/120'
[0pi/2]）
结果图形右：
0.1:
pi;
y=sin（x）;
plot（x,y,'
-sk'
holdon
pi/2;
y=x;
-b*'
）
x-x.A3/6'
[0,pi/2,0,2],2e-3,'
-gx'
holdonfplot（'
x-x.A3/6+x.A5/120'
-r.'
）holdoff
legend（'
sin（x）'
'
x'
x-A3/6+xA5/120'
2）
xlabel（'
ylabel（'
y'
title（'
Taylorapproximation'
图中红色点线为正弦曲线，蓝色的星线为一阶泰勒逼近，
绿色叉线为二阶
豪勒逼近，黑色正方形线为三阶泰勒逼近。
可见三阶泰勒逼近效果最好，泰勒级改越高，逼近效果尬婷。
3.2绘制飞机的降落曲线
一架飞机飞临北京国际机场上空时，其水平速度为540km/h,飞行高度为1000m。
飞机从距机场指挥塔的横向距离12000m处开始降落。
根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。
建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点P（x,y），则x表示飞机距指挥塔的距离，y表示飞机的飞行高度，降落曲线为
23
y（x）a0a1xa2xa3x
该函数满足条件：
y（0）0,y（12000）1000
（0）0,y'
（12000）0
（1）试利用y（x）满足的条件确定三次多项式中的四个系数；
（2）用所求出的三次多项式函数绘制出飞机降落曲线。
functionS=f（x1,y1,x2,y2,x3,y3,x4,y4）formatlong
a1=[1,x1,x1A2,x1A3];
a2=[1,x2,x2A2,x2A3];
a3=[0,1,2*x3,3*x3A2];
a4=[0,1,2*x4,3*x4A2];
a=[a1;
a2;
a3;
a4];
b=[y1;
y2;
y3;
y4];
s=a\b;
x=-12000:
250:
y=s（3）*x.A2-s（4）*x.A3
-d'
以上为M文件内容，在命令窗口键入
f（0,0,12000,1000,0,0,12000,0）
运行结桌：
ans为多项土索数
1.0e-004*
0.20833333333333
-0.00001157407407
4.1曾任英特尔公司董事长的摩尔先生早在1965年时，就观察到一件很有趣的现象：
集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。
因而发表论文，提出了大大有名的
年代
1959
1962
1963
1964
1965
增加倍数
3
5
6
摩尔定律（Moore'
sLaW,并预测未来这种增长仍会延续下去。
下面数据中，第二行数据为晶片上晶体数目在不同年代与1959年时数目
比较的倍数。
这些数据是推出摩尔定律的依据：
试从表中数据出发，推导线性拟合的函数表达式。
解法一
MATLAB：
码：
x=[1959,1962,1963,1964,1965];
y=[1,3,4,5,6];
p1=polyfit（x,y,1）
y1=polyval（p1,x）
plot（x,y1,'
-'
x,y,'
r*'
）,ylabel（'
运行结果：
p1=
1.0e+003*
0.0008-1.6255
y1=
0.81133.30194.13214.96235.7925
线性拟合的函数表达式：
Y=0.8302x-1.6255e+003
解法二：
x=[19591962196319641965];
xs=
-1.625528301891238
0.000830188679248
y=[1;
3;
4;
5;
6];
fori=1:
length（x）
forj=1:
A（i,j）=x（i）”1）;
[L,U]=lu（A'
*A）;
xs=U\（L\（A'
*y））
从而年代y与增加倍数x之间的关系为：
y=-1625.528301891238+0.830188679248x
11=
pi
12=3
13=
3.1333
14=
3.1416
（1）牛顿-莱布尼茨公式；
（2）梯形公式；
（3）辛卜生公式；
（4）复合梯形公式。
symsx
i1=int（4/（1+xA2）,x,0,1）
a=0;
b=1;
h=b-a;
i2=（4/（1+aA2）+4/（1+bA2））/2
i3=h/6*（4/（1+aA2）+4*4/（1+（a+b/2）A2）+4/（1+bA2））
M=100;
h=（b-a）/M;
i4=0;
（M-1）
x=a+h*k;
i4=i4+4/（1+x"
2）;
i4=h*（4/（1+aA2）+4/（1+bA2））/2+h*i4
牛顿莱布尼兹公式得到精确结果a采用挣形公式得到的结果比采用Simpson公式的精蒲I度鱼很，米用复祝榜形公式在步长取得越束裁小的状态下可以提高精度。
5.5求空间曲线L:
ysint
z2costsint
弧长公式为
L:
.x'
2（t）y'
2（t）z'
2（t）dt
xcost
6.1用欧拉公式和四阶龙格-库塔法分别求解下列初值问题;
0.9y
（i）y'
y（0）i;
x[0,1]
12x
2；
⑵y'
闩顶。
（1）<
1>
欧拉公式：
function[t,x]=Euler（fun,t0,tt,x0,N）h=（tt-t0）/N;
t=t0+[0:
N]'
*h;
x（1,:
）=x0'
N
f=feval（fun,t（k）,x（k,:
））;
f=f;
x（k+1,:
）=x（k,:
）+h*f;
以'
Euler.m'
保存
functionf=Euler_fun（t,x）f=0.9*x./（1+2*t）「以'
Euler_fun.m'
保存functionmain_Euler[t,x]=Euler（'
Euier_fun'
0,1,1,20）;
fh=dsolve（'
Dx=0.9*x/（1+2*t）'
x（0）=1'
8
ft（k）=t（k）;
fx（k）=subs（fh,ft（k））;
[t,x],4
main_Euler.m'
保存输入：
main_Euler
<
2>
四阶龙杯-库塔法
functionR=rk4（f,a,b,ya,N）
%y'
=f（x,y）
%a,b为左右端点
%Nf
升级会员 