数值分析课程设计报告.docx
《数值分析课程设计报告.docx》由会员分享,可在线阅读,更多相关《数值分析课程设计报告.docx(24页珍藏版)》请在冰点文库上搜索。
![数值分析课程设计报告.docx](https://file1.bingdoc.com/fileroot1/2023-5/6/64dc1e37-e72f-4d73-a579-d00459c1fd77/64dc1e37-e72f-4d73-a579-d00459c1fd771.gif)
数值分析课程设计报告
本科课程设计(报告)
数值分析课程设计
[设计题一]
设计实验验证Hilbert矩阵的病态性。
(提示:
先取x(比如x=(1,1…..1))计算出b=Hx,然后通过列主元GAUSS消去法求解Hx=b,得到近似解
比较之,并研究随着n增大,解的误差变化情况,得出结论)
设计思路:
(1)取定初值X=(1,1,1,…),计算出b=Hx,然后通过列主元GAUSS消去法求解Hx=b,得到近似解B=(b1,b2,b3…);
(2)再把近似解B,代入方程b=Hx,,解出X1,比较X与X1,判断H矩阵的病态性;
(3)比较各个H矩阵的病态性,研究它们的误差变化情况。
算法步骤:
(1)取一个三阶H矩阵,按照设计思路的方法,判断出三阶H矩阵的病态性;
(2)取一个四、五阶H矩阵,同理,判断出四阶H矩阵的病态性。
程序清单:
Xa=[1;1;1];Ha=[1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5];
Ba=Ha*Xa
Ba=[1.83;1.08;0.78];
Xa=Ha\Ba
Xb=[1;1;1;1];Hb=[1,1/2,1/3,1/4;1/2,1/3,1/4,1/5;1/3,1/4,1/5,1/6;1/4,1/5,1/6,1/7;];
Bb=Hb*Xb
Bb=[2.08;1.28;0.95;0.76];
Xb=Hb\Bb
Xc=[1;1;1;1;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)&(k>3)
disp('祝贺您!
此迭代序列收敛,且收敛速度较快')
return;
end
p=[(i-1)pchaxpchaxk]';
牛顿迭代的M文件:
function[k,xk,yk,pcha,xpcha]=newton(x0,tx,fx,n)
%k为迭代初始值,xpcha,pcha是误差范围,迭代终止条件,n控制迭代次数
x
(1)=x0;
fori=1:
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);i=i+1;
xk=x(i);yk=fun(x(i));[(i-1)xkykpchaxpcha]
if(abs(yk)k=i-1;xk=x(i);[(i-1)xkykpchaxpcha]%循环次数减一次
return;
end
end
ifi>n
k=i-1;xk=x(i);[(i-1)xkykpchaxpcha]
return;
end
[(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,而且与题目所给的
,十分接近。
牛顿迭代法:
运行结果
k=
5
xk=
1.36880810782137
yk=
0
pcha=
1.7863568394251e-015
xpcha=
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)
y=8332.5*sqrt(1-(1522.5/8332.5)^2*sin(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;%初次调用梯形公式
fori=1:
m%复化梯形公式主体
h=h/2;n=2*n;s=0;
fork=1:
n/2
x=a+h*(2*k-1);s=s+feval(fun4,x);
end
T(i+1)=T(i)/2+h*s;%最后一次调用复化梯形公式
end
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;
end
if(nargin<3)
warning(['usage:
',mfilename,'(fun,a,b[,epsilon])']);
end
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;
end
n=n+n;
h=h/2;
y2=y2+y3;
y0=y;
yiter=[yiter,y0];
end
运行过程与输出结果截图:
在命令窗口中输入:
Q1=rctrap(@fun4,0,pi/2,15)
vpa(Q1,7)
[y,n,yiter]=SimpsonInt(@fun4,0,pi/2,0.01)
vpa(y,7)
可得以下截图
结果分析:
复化梯形公式
vpa(Q1,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]
复化辛普森公式
vpa(y,7)
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是右端点值
x
(1)=x0;y=y0;%赋初值
fori=1:
n
x(i+1)=x(i)+h;
y(:
i+1)=y(:
i)+h*feval(fun,x(i),y(:
i));%欧拉公式
end
plot(x,y,'-')%作图
xlabel('x'),ylabel('y'),
(3)改进的欧拉公式:
function[x,y]=Impeuler(fun,x0,xt,y0,n,h);%定义改进的欧拉函数
x
(1)=x0;y=y0;%赋值
fori=1:
n
x(i+1)=x(i)+h;%改进欧拉公式
y1=y(:
i)+h*feval(fun,x(i),y(:
i));
y2=y(:
i)+h*feval(fun,x(i+1),y1);
y(:
i+1)=(y1+y2)/2;
end
plot(x,y,'-')%作图
xlabel('x'),ylabel('y')
(4)四阶经典龙格—库塔法:
function[x,y]=RK4(fun,y0,h,a,n);
x
(1)=a;%赋值
y(:
1)=y0;
fori=1:
n%四阶经典龙格—库塔法公式
x(i+1)=x(i)+h;
k1=feval(fun,x(i),y(:
i));
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+1)=y(:
i)+h*(k1+2*k2+2*k3+k4)/6;
end
plot(x,y,'-')%作图
xlabel('x'),ylabel('y'),
运行过程与输出结果截图:
%在命令窗口中输入,可得以下结果
[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;0],40,0.025)
[x,y]=euler(@dy,0,1,[0.5236;0],10,0.1)
[x,y]=Impeuler(@dy,0,1,[0.1745;0],20,0.05)
[x,y]=Impeuler(@dy,0,1,[0.1745;0],10,0.1)
[x,y]=Impeuler(@dy,0,1,[0.5236;0],20,0.05)
[x,y]=Impeuler(@dy,0,1,[0.5236;0],10,0.1)
[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的操作程度和对数值分析中的知识都有了极大的提高,这些或许就是我在这次课程设计中的最大收获吧。