数值分析课程设计报告Word格式文档下载.docx
《数值分析课程设计报告Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计报告Word格式文档下载.docx（24页珍藏版）》请在冰点文库上搜索。
Xc=[1;
Hc=[1,1/2,1/3,1/4,1/5;
1/2,1/3,1/4,1/5,1/6;
1/3,1/4,1/5,1/6,1/7;
1/4,1/5,1/6,1/7,1/8;
1/5,1/6,1/7,1/8,1/9;
Bc=Hc*Xc
Bc=[2.28;
1.45;
1.09;
0.89;
0.75];
Xc=Hc\Bc
运行过程与输出结果截图：
结果分析：
代入近似解Ba,Bb,Bc而得到的X如下：
Xa=
0.99000000000000
1.08000000000001
0.89999999999999
Xb=
1.27999999999997
-1.79999999999961
7.199********903
-2.79999999999936
Xc=
1.0e+002*
-0.06999999999975
1.48199999999542
-6.25799999998035
9.37999999997045
-4.53599999998560
由此分析比较初始的Xa,Xb,Xc可以看出，Hilbert矩阵是一个病态矩阵，且随着n的增大，其误差变化加大。
即，Hilbert矩阵是一个随n增大而误差偏离更大的病态矩阵。
[设计题二]
1225年，达。
芬奇研究了方程
并得到它的一个近似根
。
没有人知道他用什么方法得到它。
设计两种方法去计算，并比较这两种方法。
设计思路与算法步骤：
不动点迭代法：
（1）先把方程变换为x=20/（x.^2+2*x+10），并编译函数M文件
（2）利用不动点迭代法，确定迭代次数，进行运算，不妨设迭代次数为20，由题目可设迭代初始值X0=1
牛顿迭代法：
（1）用牛顿迭代法，得到迭代方程x1=x0－f（x0）/fˊ（x0），编译函数M文件
（2）取适当的ε，δ，和迭代次数n,进行迭代，直到符合条件时，停止运算。
（3）确定迭代初始值x0=1，设精确为0.000000001.
1、函数文件：
原函数;
functiony=fun（x）
y=x.^3+2*x.^2+10*x-20;
迭代函数：
functiony=fun1（x）
y=20/（x.^2+2*x+10）;
导函数：
functiony=dfun（x）
y=3*x.^2+4*x+10;
2、不动点迭代的M文件：
function[k,pcha,xpcha,xk]=diedai（x0,k）
%输入的量--x0是初始值,k是迭代次数
x
（1）=x0;
fori=1:
k
x（i+1）=fun1（x（i））;
%程序中调用的fun1.m为函数y=φ（x）
pcha=abs（x（i+1）-x（i））;
xpcha=pcha/（abs（x（i+1））+eps）;
%偏差计算
i=i+1;
xk=x（i）;
[（i-1）pchaxpchaxk]
end
if（pcha>
1）&
（xpcha>
0.5）&
（k>
3）%迭代收敛性判断
disp（'
请用户注意：
此迭代序列发散，请重新输入新的迭代公式'
）
return;
end
if（pcha<
0.001）&
（xpcha<
0.0000005）&
3）
祝贺您！
此迭代序列收敛，且收敛速度较快'
p=[（i-1）pchaxpchaxk]'
;
牛顿迭代的M文件：
function[k,xk,yk,pcha,xpcha]=newton（x0,tx,fx,n）
%k为迭代初始值，xpcha，pcha是误差范围，迭代终止条件，n控制迭代次数
n
x（i+1）=x（i）-fun（x（i））/（dfun（x（i））+eps）;
%牛顿迭代公式与偏差计算表达式
pcha=abs（x（i+1）-x（i））;
xpcha=pcha/（abs（x（i+1））+eps）;
yk=fun（x（i））;
[（i-1）xkykpchaxpcha]
if（abs（yk）<
fx）&
（（pcha<
tx）|（xpcha<
tx））%当偏差且自变量x的偏差小于初始精度时，算法继续。
k=i-1;
xk=x（i）;
[（i-1）xkykpchaxpcha]%循环次数减一次
end
ifi>
n
[（i-1）xkykpchaxpcha]
[（i-1）,xk,yk,pcha,xpcha]'
%在命令窗口中输入以下语句，可得以下结果:
[k,pcha,xpcha,xk]=diedai（1,20）
[k,xk,yk,pcha,xpcha]=newton（1,0.000000001,0.000000001,1000）
不动点迭代
牛顿迭代
迭代法：
运行结果
k=
20
pcha=
1.077823423845103e-007
xpcha=
7.874174939314717e-008
xk=
1.36880807468942
容易看出，进行20次迭代之后，得到的结果x=1.36880807468942，而且与题目所给的
，十分接近。
5
1.36880810782137
yk=
0
1.7863568394251e-015
1.2987398576591e-015
可知，利用牛顿迭代法进行五次迭代后，即可得到满足精度的结果，与题目给的答案接近。
总而言之，不动点迭代法与牛顿迭代法这两种方法，可见，牛顿迭代法，运算次数少，且误差小于不动点迭代的运算误差，故牛顿迭代法比较优越。
[设计题四]
地球卫星飞行的轨道是一个椭圆，椭圆的周长计算公式为
其中a为椭圆长半轴，c是椭圆中心与地球中心的距离。
令H,h分别为远地点，和近地点距离。
R=6371（km）为地球半径。
则a=（2R+H+h）/2,c=（H-h）/2.我国第一颗人造卫星H=3484km,h=439km.
（a=（2*6371+3484+439）/2=8332.5,c=（3484-439）/2）=1522.5
试通过复化梯形公式，复化辛普森公式计算卫星轨道的周长，且研究并比较这两种的方法的收敛性。
（1）、编译被积函数M文件，被积函数表达式为：
y=8332.5*sqrt（1-（1522.5/8332.5）^2*sin（x））;
（2）、根据复化梯形公式，编译算法。
不妨设递归次数为15次。
（由于复化梯形公式，的两个端点值，所使用的公式不同，故要分开编写循环体）
（3）、根据复化辛普森公式，
，编译算法，不妨设递归次数为15次。
被积函数：
functiony=fun4（x）
复化梯形公式;
functionT=rctrap（fun4,a,b,m）%fun4表示被积函数，a，b是被积函数上下限，m是递归次数
n=1;
h=b-a;
T=zeros（1,m+1）;
x=a;
%赋值计算
T
（1）=h*（feval（fun4,a）+feval（fun4,b））/2;
%初次调用梯形公式
m%复化梯形公式主体
h=h/2;
n=2*n;
s=0;
fork=1:
n/2
x=a+h*（2*k-1）;
s=s+feval（fun4,x）;
T（i+1）=T（i）/2+h*s;
%最后一次调用复化梯形公式
T=T（1:
m）;
复化辛普森公式：
function[y,n,yiter]=SimpsonInt（fun,a,b,epsilon）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Function:
Simpson积分
%%Argument:
fun是被积函数；
a,b是积分上下限；
epsilon是误差容限
%%Return:
y为积分结果；
将[a,b]平均n等分；
yiter迭代的积分值
if（nargin<
4）
epsilon=1e-6;
warning（['
usage:
'
mfilename,'
（fun,a,b[,epsilon]）'
]）;
n=2;
h=（b-a）/4;
y1=feval（fun,a）+feval（fun,b）;
y2=feval（fun,（a+b）/2）;
y0=（b-a）/6*（y1+4*y2）;
yiter=y0;
while1
y3=sum（feval（fun,a+（1:
2:
2*n-1）*h））;
y=h/3*（y1+2*y2+4*y3）;
ifabs（y-y0）<
15*epsilon
break;
n=n+n;
h=h/2;
y2=y2+y3;
y0=y;
yiter=[yiter,y0];
在命令窗口中输入：
Q1=rctrap（@fun4,0,pi/2,15）
vpa（Q1,7）
[y,n,yiter]=SimpsonInt（@fun4,0,pi/2,0.01）
vpa（y,7）
可得以下截图
复化梯形公式
ans=
[12978.49,12955.87,12950.43,12949.09,12948.75,12948.67,12948.65,12948.64,12948.64,12948.64,12948.64,12948.64,12948.64,12948.64,12948.64]
复化辛普森公式
ans=12948.64
可见，利用复化梯形公式计算得到的周长为12948.64,而利用复化辛普森公式计算得到的周长为12948.64时，仅递归了4次。
由此可见，复化辛普森公式的收敛速度比复化梯形公式快。
[设计题五]
给定单摆方程初值问题
其中g=9.8,l=25.
取初始偏离角度
其精确解为
分别对上述两种情况按照下列方法求出其数值解，比较各方法的优缺点，并将计算结果与精确解做比较（列表、画图）。
（方案I）欧拉法，步长h=0.025,h=0.1；
（方案II）改进的欧拉法，步长h=0.05,h=0.1；
（方案III）四阶经典龙格—库塔法，步长h=0.1。
对于高阶常微分方程的解法：
1、把二阶微分方程化成标准型dt/dx=-g/l*sinxt（0）=0x（0）=θx=θ;
t=x’
2、把函数编译成矩阵表达式
3、解出精确解
4、利用欧拉法，改进的欧拉法，四阶龙格库塔法，求解
5、在相同的步长内，比较三种方法和精确解的误差大小
函数M文件:
functionf=dy（x,y）
f=[y
（2）;
-9.8*sin（y
（1））/25];
（1）精确解：
y=0.1745*cos（9.8*x/25）
plot（x,y,'
-'
）;
xlabel（'
x'
）,ylabel（'
y'
），
（2）欧拉法：
function[x,y]=Euler（fun,x0,xt,y0,n,h）;
%xt是右端点值
y=y0;
%赋初值
x（i+1）=x（i）+h;
y（:
i+1）=y（:
i）+h*feval（fun,x（i）,y（:
i））;
%欧拉公式
）%作图
）,
（3）改进的欧拉公式：
function[x,y]=Impeuler（fun,x0,xt,y0,n,h）;
%定义改进的欧拉函数
%赋值
n
%改进欧拉公式
y1=y（:
y2=y（:
i）+h*feval（fun,x（i+1）,y1）;
i+1）=（y1+y2）/2;
）
（4）四阶经典龙格—库塔法：
function[x,y]=RK4（fun,y0,h,a,n）;
x
（1）=a;
%赋值
1）=y0;
n%四阶经典龙格—库塔法公式
x（i+1）=x（i）+h;
k1=feval（fun,x（i）,y（:
k2=feval（fun,x（i）+h/2,y（:
i）+h*k1/2）;
k3=feval（fun,x（i）+h/2,y（:
i）+h*k2/2）;
k4=feval（fun,x（i）+h,y（:
i）+h*k2）;
y（:
i）+h*（k1+2*k2+2*k3+k4）/6;
%在命令窗口中输入，可得以下结果
[x,y]=euler（@dy,0,1,[0.1745;
0],10,0.1）
[x,y]=euler（@dy,0,1,[0.1745;
0],40,0.025）
[x,y]=euler（@dy,0,1,[0.5236;
[x,y]=Impeuler（@dy,0,1,[0.1745;
0],20,0.05）
[x,y]=Impeuler（@dy,0,1,[0.1745;
[x,y]=Impeuler（@dy,0,1,[0.5236;
[x,y]=RK4（@dy,[0.5236;
0],0.1,0,10）
[x,y]=RK4（@dy,[0.1745;
0],0.1,0,10）
y=0.1745*cos（sqrt（9.8/25）*x）
y=0.5236*cos（sqrt（9.8/25）*x）
截图：
1、精确解：
（方案I）欧拉法：
（1）步长h=0.025，
（2）h=0.1
（3）h=0.025
（4）h=0.1
（方案II）改进的欧拉法：
（1）步长h=0.05,
（2）h=0.1,
（1）h=0.1,
（2）h=0.05,
（方案III）四阶经典龙格—库塔法:
（1）
（2）
1、当
，h=0.1时，
从上诉的程序运行过程和结果可以看出，利用改进的欧拉法和经典四阶龙格库塔公式解出的值，比较接近精确解，而欧拉法比较远离精确解。
且经典四阶龙格库塔公式收敛相对快一点。
2、当
，h=0.1时，
得出的结论，与
时，相同。
总而言之，利用四阶龙格库塔方法所得的结果比较具有优越性，且收敛快。
而欧拉算法，比较简单，运算比较快，改进的欧拉法处于两者之间
心得体会
这次课程设计的时间不算长，可是我从中学习到了很多东西。
由于有选择数学软件这个课程，所以对于matlab这个软件的认识还是有一定的基础。
不过，由于实际操作比较小，所以处理起这个课程设计还是比较吃力。
从一开始，我对于matlab的一无所知，到现在可以利用这个软件对一些常用的算法，加于研究和设计，途中经历了很多困难。
比如，如何使用for语句，矩阵的运算操作等等，都不太清楚。
具体的，例如设计题二。
虽然利用数值分析上所学到的知识，可以从理论上解决这个问题。
但是要利用matlab来处理的话，我顿时觉得无从下手，因为，要设计一个算法，毕竟涉及很多细小的问题，比如对于不动点迭代，需要如何设计变量，需要什么变量，for循环中的fun改怎么使用，怎么改变下标。
迭代过程的终止条件改如何设计等等。
在遇到这些困难的时候，我要做的最多的就是，不断的搜索和阅读材料。
从实例中总结和学习他人的方法，然后挪用到自己的程序算法中，不断的调试，修改，最后才把自己的算法设计出来。
在完成那一刻，心中便产生起一种莫名的成就感，也许这些就是我们在成功解决问题后，所得到的心里上最大的回报吧。
虽然，当你对一个无知的事物加以学习的过程是比较繁琐的。
但是，当你弄懂，和能够分析它的时候，你会发现，你能感受到一种发自内心无的成就感冲击着你的大脑。
只有尝试过，你才知道这个过程的乐趣。
通过这个程序设计，我对于matlab软件的使用，无疑更加熟练了。
而且对于数值分析这个课程，我也有了更深刻的体会。
通过这个课程设计，我认识到，为什么说，计算方法是为各种数学问题的数值解答研究提供最有效的算法。
因为利用计算方法中的知识，我们可以再误差允许的范围内，用最快，最准确的方法去解决我们生活中遇到的问题。
例如，卫星运行轨迹的周长运算，虽然可以得到精确的运算表达式，但是按部就班的去积分运算的话，无疑是一件复杂的工程。
而利用数值分析中的辛普森和梯形公式去计算，虽然有误差，但是误差范围可以控制，而且计算过程快捷。
这就是数值分析在实际生活中存在的意义吧。
总而言子，通过这次课程设计，我深刻认识到，我今天认真的进行课程设计，学会脚踏实地迈开这一步，就是为明天能稳健地在社会大潮中奔跑打下坚实的基础。
在这次课程设计中，我的处理问题能力，和信息检索能力，对matlab的操作程度和对数值分析中的知识都有了极大的提高，这些或许就是我在这次课程设计中的最大收获吧。
升级会员 