数值分析课程设计.docx
《数值分析课程设计.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计.docx（25页珍藏版）》请在冰豆网上搜索。
数值分析课程设计
数值分析
课程设计报告
专业：
班级：
学号：
姓名：
1．1水手、猴子和椰子问题：
五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。
由于旅途的颠簸，大家都很疲惫，很快就入睡了。
第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？
试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题（15621）。
分析：
定义数组x（）用于存储每天的椰子数
通过给最后一天的椰子数x（7）附初值根据公式“x（n）=5*x（n+1）/4+1”逆向递推前一天的椰子数
如果推出某一天的椰子数不为整数（即：
mod（5*x（n+1）,4）~=0），说明初值x（7）不正确，则进行x（7）=x（7）+1运算来改变初值，直到符合要求为止。
最终一次性输x（），6次分配的前后值。
程序代码：
functionyezi1（m）
x（7）=0;
whilex（7）x（7）=m;
x（6）=5*x（7）+1;
n=5;
while（mod（5*x（n+1）,4）==0）&（n>0）
x（n）=5*x（n+1）/4+1;
n=n-1;
end
ifn~=0
m=x（7）+1;
end
end
disp（x）
指令窗口输入：
>>m=1;yezi1（m）
回车输出：
156211249699967996639651161023
本程序，首先赋值m=1（m为最后分完时每堆椰子数量），程序自身会测试m的可行性，如不满足要求，程序自动增大m值，直到出现第一个结果，此时m=1023，总椰子数为15261.
如果输入m初值>1023,程序会自动找到下一个符合值2047，…………
2．1用高斯消元法的消元过程作矩阵分解。
设
消元过程可将矩阵A化为上三角矩阵U，试求出消元过程所用的乘数
、
、
并以如下格式构造下三角矩阵L和上三角矩阵U
验证：
矩阵A可以分解为L和U的乘积，即A=LU。
分析：
首先用矩阵A的各阶主子式值判断矩阵能否进行LU分解，在能够分解的前提下，使用高斯消去法，将将消去时的两行间的倍数m附入L（m=A（k,p）/A（p,p）;L（p,p）=1;L（k,p）=m;），再将消去后的A附入U，即得L,U。
程序代码：
function[L,U]gaussLU（A）
[nn]=size（A）;RA=rank（A）;
ifRA~=n
disp（'因为A的n阶行列式hl等于零，所以A不能进行LU分解。
'）;
hl=det（A）;
return
end
ifRA==n
forp=1:
n
h（p）=det（A（1:
p,1:
p））;
end
hl=h（1:
n）;
fori=1:
n
ifh（1,i）==0
disp（'因为A的r阶主子式等于零，所以A不能进行LU分解。
'）
return
end
end
ifh（1,i）~=0
disp（'因为A的各阶主子式都不等于零，所以A能进行LU分解。
L与U的值如下：
'）
L=zeros（n,n）;L（n,n）=1;
U=zeros（n,n）;
forp=1:
n-1
fork=p+1:
n
m=A（k,p）/A（p,p）;
L（p,p）=1;L（k,p）=m;
A（k,p:
n）=A（k,p:
n）-m*A（p,p:
n）;
end
end
U=A（1:
n,1:
n）;
end
end
指令窗口输入：
A=[2023;181;235];[L,U]=gaussLU（A）
程序输出：
因为A的各阶主子式都不等于零，所以A能进行LU分解。
L与U的值如下：
L=
1.000000
0.05001.00000
0.10000.35441.0000
U=
20.00002.00003.0000
07.90000.8500
004.3987
2．2用矩阵分解方法求上题中A的逆矩阵。
记
分别求解方程组
由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。
对第一个方程组，由于A=LU，所以先求解下三角方程组
再求解上三角方程组
，则可得逆矩阵的第一列列向量；类似可解第二、第三方程组，得逆矩阵的第二列列向量的第三列列向量。
由三个列向量拼装可得逆矩阵
。
分析：
我觉得题目中已经分析的很全面了，此题没有必要再多此一举了。
程序代码：
A=[20,2,3;1,8,1;2,-3,15];
[L,U]=gaussLU（A）;
b1=[1;0;0];b2=[0;1;0];b3=[0;0;1];
Y1=L\b1;
Y2=L\b2;
Y3=L\b3;
X1=U\Y1;
X2=U\Y2;
X3=U\Y3;
X=[X1,X2,X3]%X
输出结果：
X=
0.0517-0.0164-0.0093
-0.00550.1237-0.0072
-0.00800.02690.0665
实验三
3．1用泰勒级数的有限项逼近正弦函数
用计算机绘出上面四个函数的图形。
分析：
直接设定变量范围、函数，然后作图即可，注意图像可能会重叠，最好用不同线条，区别大一点的颜色绘制。
程序代码：
x=0:
pi/100:
pi;
x1=0:
pi/100:
pi/2;
y0=sin（x）;
y1=x1;
y2=x1-x1.^3/6;
y3=x1-x1.^3/6+x1.^5/120;
plot（x,y0,':
k',x1,y1,'-g',x1,y2,'-b',x1,y3,'-r'）
输出图像：
3．2绘制飞机的降落曲线
一架飞机飞临北京国际机场上空时，其水平速度为540km/h，飞行高度为1000m。
飞机从距机场指挥塔的横向距离12000m处开始降落。
根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。
建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点
，则
表示飞机距指挥塔的距离，
表示飞机的飞行高度，降落曲线为
该函数满足条件：
（1）试利用
满足的条件确定三次多项式中的四个系数；
（2）用所求出的三次多项式函数绘制出飞机降落曲线。
分析：
通过y（x）满足的条件带入降落曲线方程，即可得到一个四元一次方程组。
将方程组视为矩阵。
由于前面的题目用过了高斯消去法，本体继续沿用此方法，通过高斯消去法，求解方程组，得到a0,a1,a2,a3的值，带回y中，运用plot函数，输出x，y的关系曲线。
程序代码：
functiona=jiangl（x1,x2,x3,x4,b）
A=[1x1x1^2x1^3;1x2x2^2x2^3;012*x33*x3*x3;012*x43*x4*x4];
B=[Ab];
a=zeros（4,1）;C=zeros（1,5）;
forp=1:
3
[Y,j]=max（abs（B（p:
4,p）））;C=B（p,:
）;
B（p,:
）=B（j+p-1,:
）;B（j+p-1,:
）=C;
fork=p+1:
4
m=B（k,p）/B（p,p）;
B（k,p:
5）=B（k,p:
5）-m*B（p,p:
5）;
end
end
b=B（1:
4,5）;A=B（1:
4,1:
4）;a（4）=b（4）/A（4,4）;
forq=3:
-1:
1
a（q）=（b（q）-sum（A（q,q+1:
4）*a（q+1:
4）））/A（q,q）;
end
在指令窗口输入：
>>x1=0;x2=12000;x3=0;x4=12000;b=[0;1000;0;0];a=jiangl（x1,x2,x3,x4,b）
回车输出：
a=
1.0e-004*
0
0
0.2083
-0.0000
>>x=0:
50:
12000;y=a
（1）+x*a
（2）+x.^2*a（3）+x.^3*a（4）;plot（x,y）
回车，输出图像：
4．1曾任英特尔公司董事长的摩尔先生早在1965年时，就观察到一件很有趣的现象：
集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。
因而发表论文，提出了大大有名的摩尔定律（Moore’sLaw），并预测未来这种增长仍会延续下去。
下面数据中，第二行数据为晶片上晶体数目在不同年代与1959年时数目比较的倍数。
这些数据是推出摩尔定律的依据：
年代
1959
1962
1963
1964
1965
增加倍数
1
3
4
5
6
试从表中数据出发，推导线性拟合的函数表达式。
分析：
通过表中数据绘图，可发现趋向于一次线性方程，正好本题要求推到线性拟合函数，直接使用一次的最小二乘法。
程序代码：
functionzxec（x,y）
y=y';
fori=1:
length（x）
forj=1:
2
A（i,j）=x（i）^（j-1）;
end
end
[L,U]=lu（A'*A）;
alpha=U\（L\（A'*y））
在指令窗口输入：
x=[1959,1962,1963,1964,1965];y=[1,3,4,5,6];zxec（x,y）
程序输出：
alpha=
1.0e+003*
-1.6255
0.0008
其中alpha
（1）=-1625.5是函数交y轴的焦点值，0.8是一阶函数斜率
5．1用几种不同的方法求积分
的值。
（1）牛顿-莱布尼茨公式；
（2）梯形公式；（3）辛卜生公式；（4）复合梯形公式。
程序代码：
（1）牛顿-莱布尼茨
由于
=4arctgx|0到1
输入内容：
>>4*atan
（1）
程序输出：
ans=
3.1416
（2）（3）（4）
定义
functions=f（x）
s=4/（1+x^2）;
输入：
a=0;b=1*
（2）梯形公式
接着*式输入：
（1/2）*（b-a）*（f（b）+f（a））
输出结果：
ans=
3
（3）辛卜生公式
接着*式输入：
（1/6）*（b-a）*（f（a）+4*f（（a+b）/2）+f（b））
输出结果：
ans=
3.1333
（4）复合梯形公式
接着*式输入：
T=0;
forn=0:
0.1:
0.9
m=0.1/2*（f（n）+f（n+0.1））;T=T+m;
end
T
输出结果：
T=
3.1399
6．1用欧拉公式和四阶龙格-库塔法分别求解下列初值问题；
分析：
此题只需建立欧拉和龙格库塔函数，带入求解即可。
由于题目没有给出步长，我们取步长h=0.1
程序代码：
（1）建立名为f.m的文件：
functionr=f（t,y）
r=（0.9*y）/（1+2*t）;
Euler方法
functioneuler（a,b,h,alpha）
t=a:
h:
b;
n=（b-a）/h;
y
（1）=alpha;
fori=2:
n+1
y（i）=y（i-1）+h*f（t（i-1）,y（i-1））;
sprintf（'t=%3.1f,y=%9.7f',t（i）,y（i））
end
在命令窗口输入：
euler（0,1,0.1,1）
回车得到：
ans=
t=0.1,y=1.0900000
t=0.2,y=1.1717500
t=0.3,y=1.2470768
t=0.4,y=1.3172249
t=0.5,y=1.3830861
t=0.6,y=1.4453250
t=0.7,y=1.5044519
t=0.8,y=1.5608688
t=0.9,y=1.6148989
t=1.0,y=1.6668064
四阶Runge-Kutta方法
functionrk4（a,b,h,alpha）
t=a:
h:
b;
n=（b-a）/h;
y
（1）=alpha;
fori=2:
n+1
k1=f（t（i-1）,y（i-1））;
k2=f（t（i-1）+h/2,y（i-1）+h/2*k1）;
k3=f（t（i-1）+h/2,y（i-1）+h/2*k2）;
k4=f（t（i-1）+h,y（i-1）+h*k3）;
y（i）=y（i-1）+h/6*（k1+2*k2+2*k3+k4）;
sprintf（'t=%3.1f,y=%9.7f',t（i）,y（i））
end
在命令窗口输入：
rk4（0,1,0.1,1）
回车得到：
ans=
t=0.1,y=1.0855051
t=0.2,y=1.1634777
t=0.3,y=1.2355334
t=0.4,y=1.3027862
t=0.5,y=1.3660420
t=0.6,y=1.4259056
t=0.7,y=1.4828446
t=0.8,y=1.5372290
t=0.9,y=1.5893579
t=1.0,y=1.6394763
同理：
（2）建立名为f.m的文件：
functionr=f（t,y）
r=-（t*y）/（1+t^2）;
在命令窗口分别输入：
euler（0,1,0.1,2）
回车
rk4（0,1,0.1,2）
回车
结果如下表
Euler方法
四阶Runge-Kutta方法
t=0.1,y=2.0000000
t=0.2,y=1.9801980
t=0.3,y=1.9421173
t=0.4,y=1.8886645
t=0.5,y=1.8235382
t=0.6,y=1.7505966
t=0.7,y=1.6733644
t=0.8,y=1.5947500
t=0.9,y=1.5169573
t=1.0,y=1.4415285
t=0.1,y=1.9900743
t=0.2,y=1.9611612
t=0.3,y=1.9156523
t=0.4,y=1.8569530
t=0.5,y=1.7888539
t=0.6,y=1.7149854
t=0.7,y=1.6384634
t=0.8,y=1.5617372
t=0.9,y=1.4865879
t=1.0,y=1.4142132
7.1用二分法求下列方程在指定区间[a，b]上的实根近似值：
（要求误差不超过0.01）
（1）x-sinx-1=0，[a,b]=[1,3]；
（2）xsinx=1，[a,b]=[0,1.5].
分析：
直接将方程带入二分法函数中求解即可。
程序代码：
functionx=eff（fun,a,b,eps）
ifnargin<4
eps=1e-5;
end
fa=feval（fun,a）;fb=feval（fun,b）;
iffa*fb>0
disp（'[a,b]不是有根区间，请重新调整'）;
return;
end
whileabs（b-a）/2>eps
x=（a+b）/2;fx=feval（fun,x）;
iffx*fa<0
b=x;fb=fx;
else
a=x;fa=fx;
end
end
x=（a+b）/2;
在命令窗口输入：
（1）fun=inline（'x-sin（x）-1'）;x=eff（fun,1,3,0.01）
输出结果
x=
1.9297
（2）fun=inline（'x*sin（x）-1'）;x=eff（fun,0,1.5,0.01）
输出结果：
x=
1.1191
8.1分别用直接法、雅可比迭代法、赛德尔迭代法求解线性方程组AX=b。
（1）
，
（2）
，
分析：
首先通过diag函数给矩阵A附上初值，然后分别调用雅可比和赛德尔函数即可
程序：
（1）
M=ones（9,9）;
N=4*ones（10,10）;
A=diag（diag（N））+diag（diag（M）,1）+diag（diag（M）,-1）;
b=[2;1;1;1;1;1;1;1;1;2];
1.用直接法输入：
formatlong
A\b
程序输出：
ans=
0.479271991911021
0.082912032355915
0.189********5319
0.160768452982811
0.167846309403438
0.167846309403438
0.160768452982811
0.189********5319
0.082912032355915
0.479271991911021
2.用Jacobi法
程序代码：
functionJacobi（A,b,X0,P,error,maxl）
[nn]=size（A）;
X=zeros（n,1）;
fork=1:
maxl
forj=1:
n
X（j）=（b（j）-A（j,[1:
j-1,j+1:
n]）*X0（[1:
j-1,j+1:
n]））/A（j,j）;
end
errX=norm（X-X0,P）;
X0=X;
if（errXdisp（'迭代次数k，近似解X分别是：
'）
k
X
return
end
end
if（errX>=error）
disp（'请注意：
Jacobi迭代次数已经超过最大迭代次数maxl.'）
end
在命令窗口输入：
X0=[0;0;0;0;0;0;0;0;0;0];
Jacobi（A,b,X0,inf,0.001,100）;
程序输出：
迭代次数k，精确解X1和近似解X分别是：
k=
9
X=
0.47933959960938
.0831********
0.189********852
0.16108322143555
0.16813659667969
0.16813659667969
0.16108322143555
0.189********852
.0831********
0.47933959960938
用赛德尔迭代法：
程序代码：
functionGS（A,b,X0,P,error,maxl）
[nn]=size（A）;
X=zeros（n,1）;
fork=1:
maxl
forj=1:
n
XX=0;
fori=1:
n
ifiXX=XX+A（j,i）*X（i）;
end
ifi>j
XX=XX+A（j,i）*X0（i）;
end
end
X（j）=（b（j）-XX）/A（j,j）;
end
errX=norm（X-X0,P）;
X0=X;
if（errXdisp（'迭代次数k，近似解X分别是：
'）
k
X
return
end
end
if（errX>=error）
disp（'请注意：
Gauss-Seidel迭代次数已经超过最大迭代次数maxl.'）
end
在命令窗口输入：
X0=[0;0;0;0;0;0;0;0;0;0];
GS（A,b,X0,inf,0.001,100）;
输出结果：
k=
6
X=
0.47925853729248
0.08289575576782
0.189********913
0.16064277291298
0.16802297346294
0.16770140314475
0.16085489129182
0.189********114
0.08292683955369
0.47926829011158
（2）
U=ones（18,18）;
V=-2*ones（19,19）;
W=5*ones（20,20）;
A=diag（diag（W））+diag（diag（V）,1）+diag（diag（U）,2）+diag（diag（V）,-1）+diag（diag（U）,-2）;
b=[1;zeros（19,1）];
X0=zeros（20,1）;
同上获得如下解：
直接法
雅可比
赛德尔
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
0.242121373548085
.0943********
.021*********
.0313********
-0.006942557862112
0.004886927715767
0.003586117214438
0.000214823135181
-0.000784526508185
-0.000357861979163
0.000050437638520
0.000106309021306
0.000029232399870
-0.000014343050585
-0.000012667696430
-0.000001464361499
0.000002491258113
0.000001311981447
-0.000000093354238
-0.000000299737984
1.0e+004*-0.28210827857172
1.0e+004*0.48204640630127
1.0e+004*-0.69195375595679
1.0e+004*0.88126694207517
1.0e+004*-1.05291766695873
1.0e+004*1.20130020954019
1.0e+004*-1.32375170735604
1.0e+004*1.41753761486760
1.0e+004*-1.48072923013544
1.0e+004*1.51202411408633
1.0e+004*-1.51082661747071
1.0e+004*1.47723879582980
1.0e+004*-1.41205236153589
1.0e+004*1.31674011324409
1.0e+004*-1.19335996640489
1.0e+004*1.04472588730159
1.0e+004*-0.87353157471502
1.0e+004*0.68530915502268
1.0e+004*-0.47711982856112
1.0e+004*0.27913406015473
0.24229171609600
0.09477250842624
.021*********
.0314********
-0.0071474575114
升级会员 