数值分析课程设计.docx
《数值分析课程设计.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计.docx（30页珍藏版）》请在冰豆网上搜索。
数值分析课程设计
数值分析课程设计报告
设计题4、6、8、11、12
姓名
学号
院系数学与信息科学学院
专业信息与计算科学
年级2013级
指导教师
2016年04月25日
目录
设计题四1
龙格现象实验1
1.1问题分析与设计思路1
在计算方法中，有利用多项式对某一函数的近似逼近，这样，利用多项式就可以计算相应的函数值。
例如，在事先不知道某一函数的具体形式的情况下，只能测量得知某一些分散的函数值。
例如我们不知道气温随日期变化的具体函数关系，但是我们可以测量一些孤立的日期的气温值，并假定此气温随日期变化的函数满足某一多项式。
这样，利用已经测的数据，应用待定系数法便可以求得一个多项式函数。
应用此函数就可以计算或者说预测其他日期的气温值。
一般情况下，多项式的次数越多，需要的数据就越多，而预测也就越准确，这就是龙格现象。
1
1.2程序清单1
（1）编写拉格朗日插值函数，将其存到当前路径的M文件中：
1
functiony=lagrange（x0,y0,x）1
n=length（x0）;m=length（x）;1
fori=1:
m1
z=x（i）;1
L=0.0;1
forj=1:
n1
T=1.0;1
fork=1:
n1
ifk~=j1
T=T*（z-x0（k））/（x0（j）-x0（k））;1
end1
end1
L=T*y0（j）+L;1
end1
y（i）=L;1
end2
（2）分别取不同的n值，作出相对应n值的插值多项式的曲线图。
2
n=4时，2
x0=-5:
10/4:
5;2
y0=1./（1+x0.^2）;2
x=-5:
0.1:
5;2
y=lagrange（x0,y0,x）;2
y1=1./（1+x.^2）;2
plot（x,y,'-k'）2
holdon2
plot（x,y,'-.r'）2
n=6时，2
x0=-5:
10/6:
5;2
y0=1./（1+x0.^2）;2
x=-5:
0.1:
5;2
y=lagrange（x0,y0,x）;2
>>y1=1./（1+x.^2）;2
>>plot（x,y1,'-k'）2
>>holdon2
>>plot（x,y,'--h'）2
n=8时，2
x0=-5:
10/8:
5;2
y0=1./（1+x0.^2）;2
x=-5:
0.1:
5;2
y=lagrange（x0,y0,x）;2
y1=1./（1+x.^2）;2
plot（x,y1,'-k'）2
holdon2
plot（x,y,'--g'）2
n=10时，2
x0=-5:
1:
5;2
y0=1./（1+x0.^2）;2
x=-5:
0.1:
5;2
y=lagrange（x0,y0,x）;3
y1=1./（1+x.^2）;3
plot（x,y1,'-k'）3
holdon3
plot（x,y,'--m'）3
1.3结果分析3
1.4设计总结3
数值分析课程设计总结22
设计题四
龙格现象实验
1.1问题分析与设计思路
在计算方法中，有利用多项式对某一函数的近似逼近，这样，利用多项式就可以计算相应的函数值。
例如，在事先不知道某一函数的具体形式的情况下，只能测量得知某一些分散的函数值。
例如我们不知道气温随日期变化的具体函数关系，但是我们可以测量一些孤立的日期的气温值，并假定此气温随日期变化的函数满足某一多项式。
这样，利用已经测的数据，应用待定系数法便可以求得一个多项式函数。
应用此函数就可以计算或者说预测其他日期的气温值。
一般情况下，多项式的次数越多，需要的数据就越多，而预测也就越准确，这就是龙格现象。
对于函数
进行拉格朗日插值，取不同的节点数n，在区间[-5,5]上取等距间隔的节点为插值点，把
和插值多项式的曲线画在同一张图上进行比较。
1.2程序清单
（1）编写拉格朗日插值函数，将其存到当前路径的M文件中：
functiony=lagrange（x0,y0,x）
n=length（x0）;m=length（x）;
fori=1:
m
z=x（i）;
L=0.0;
forj=1:
n
T=1.0;
fork=1:
n
ifk~=j
T=T*（z-x0（k））/（x0（j）-x0（k））;
end
end
L=T*y0（j）+L;
end
y（i）=L;
end
（2）分别取不同的n值，作出相对应n值的插值多项式的曲线图。
n=4时，
x0=-5:
10/4:
5;
y0=1./（1+x0.^2）;
x=-5:
0.1:
5;
y=lagrange（x0,y0,x）;
y1=1./（1+x.^2）;
plot（x,y,'-k'）
holdon
plot（x,y,'-.r'）
n=6时，
x0=-5:
10/6:
5;
y0=1./（1+x0.^2）;
x=-5:
0.1:
5;
y=lagrange（x0,y0,x）;
>>y1=1./（1+x.^2）;
>>plot（x,y1,'-k'）
>>holdon
>>plot（x,y,'--h'）
n=8时，
x0=-5:
10/8:
5;
y0=1./（1+x0.^2）;
x=-5:
0.1:
5;
y=lagrange（x0,y0,x）;
y1=1./（1+x.^2）;
plot（x,y1,'-k'）
holdon
plot（x,y,'--g'）
n=10时，
x0=-5:
1:
5;
y0=1./（1+x0.^2）;
x=-5:
0.1:
5;
y=lagrange（x0,y0,x）;
y1=1./（1+x.^2）;
plot（x,y1,'-k'）
holdon
plot（x,y,'--m'）
1.3结果分析
1.4设计总结
上述现象及结果告诉我们，在进行数值计算时，并不是插值多项式的次数越高（即插值节点越多），插值效果越好，精度也不一定随次数的提高而升高，这种现象称为Runge现象。
从数值计算上可解释为高次插值多项式的计算会带来舍入误差的增大，从而引起计算失真。
因此，实际应用做插值时一般只用一次、二次最多用三次插值多项式。
设计题六
函数逼近Matlab实现及其应用
2.1问题思路与设计思路
在数值计算中经常要计算函数值，如计算机中计算基本初等函数及其他特殊函数；当函数只在有限点集上给定函数值，要在包含该点集的区间上用公式给出函数的简单表达式，这些都涉及在区间[a,b]上用简单函数逼近已知复杂函数的问题，这就是函数逼近问题。
常用的有最佳平方逼近、最小二乘法。
2.2程序清单
（1）最佳平方逼近：
legendre（N）函数：
functionp=legendre（N）
symstx;%定义符号变量tx
forn=1:
N
pp（n）=diff（（t^2-1）^（n-1）,n-1）;%diff函数，求函数的n阶导数
Q（n）=2^（n-1）*prod（[1:
n-1]）;%prod函数，计算数组元素的连乘积
end
pp
（1）=1;
Q=sym（Q）;
p=pp*（inv（diag（Q）））;
用M文件建立被逼近函数：
functionF=creat（x）
n=length（x）;
F=x.*cos（x（1:
n））;
区间变换函数程序：
functionf=convert（a,b,F）
symsxt;
s=2\（（b-a）*t+a+b）;
f=subs（F,x,s）;
变步长复化梯形求积公式：
functionI=tx（g）
m=1;
h=1-（-1）;
T=zeros（1,100）;
T
（1）=h*（feval（g,-1）+feval（g,1））/2;
i=1;
whilei<100
m=2*m;
h=h/2;
s=0;
fork=1:
m/2
x=-1+h*（2*k-1）;
s=s+feval（g,x）;
end
T（i+1）=T（i）/2+h*s;
ifabs（T（i+1）-T（i））<0.00001
I=T（i+1）;
break;
end
i=i+1;
end
最佳平法逼近函数leastp：
function[cs]=leastp（a,b,N）
symstx;
F=creat（x）;
p=legendre（N）;
f=convert（a,b,F）;
f=p*diag（f）;
fori=1:
N
g=inline（f（i））;
I=tx（g）;
u（i）=I;
Q（i）=2\（2*（i-1）+1）;
end
Q=sym（Q）;
c=double（u*diag（Q））;
S=c*p';
s=subs（S,t,（2*x-a-b）/（b-a））;
subplot（211）,
ezplot（s,[a:
0.01:
b]）;
subplot（212）,
ezplot（F,[a,b]）;
在命令窗口输入以下代码：
a=0;b=4;N=2;
[cs]=leastp（a,b,N）;
a=0;b=4;N=4;
[cs]=leastp（a,b,N）;
a=0;b=4;N=7;
[cs]=leastp（a,b,N）;
（2）最小二乘法
x=[1:
1:
30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit（x,y,3）
a2=polyfit（x,y,9）
a3=polyfit（x,y,15）
b1=polyval（a1,x）
b2=polyval（a2,x）
b3=polyval（a3,x）
r1=sum（（y-b1）.^2）
r2=sum（（y-b2）.^2）
r3=sum（（y-b3）.^2）
2.3结果分析
（1）最佳平方逼近
N=2时，
N=4时，
N=7时，
（2）最小二乘法
在命令窗口输入：
plot（x,y,'*'）
holdonplot（x,b1,'r'）
holdonplot（x,b2,'g'）
holdonplot（x,b3,'b:
o'）
运行结果：
2.4设计总结
通过本次实验，我不仅对从前所学的数值线性代数的知识有了一定的复习和理解，包括最小二乘法，拟合数据，拟合多项式等有了更好的把握，并且练习了如何更准确和熟练地使用Matlab 软件，该款软件非常适用于计算量大的问题，不仅简单易学而且编程效率高。
它包含很多软件包，方便又实用。
当然，我在实验过程中也遇到了不少问题，但是在同学的帮助下都得到了妥善的解决，我认为，此次最大的收获并不是解决了某一个问题和认识了一个软件，而是教会我们如何自己去了解和学习使用计算工具，这对我们今后的学习和工作有很大的帮助。
设计题八
常见数值积分方法的Matlab实现及其应用
3.1问题思路与设计思路
实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系。
依据人们所熟知的微积分基本定理，对于积分只要找到被积函数f（x）的原函数F（x），便有牛顿—莱布尼茨公式：
但实际使用这种求积方法往往有困难，因为有大量的被积函数其原函数不能用初等函数表达，故不能用上述公式计算，因此有必要研究积分的数值计算问题。
但实际使用这种求积方法往往有困难，因为有大量的被积函数其原函数不能用初等函数表达，故不能用上述公式计算，因此有必要研究积分的数值计算问题。
数值求积方法是近似方法，为要保证精度，我们自然希望求积公式能对尽可能多的函数准确的成立，这就提出了所谓代数精度的概念：
如果某个求积公式对于次数不超过m的多项式均能准确的成立，但对于m+1次多项式就不准确成立，则称该求积公式具有m次代数精度。
构造数值积分公式最通常的方法是用积分区间上的n次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式。
3.2程序清单
复合辛普森公式求解：
functions=xinpusen（fun,a,b,n）
h=（b-a）/n
s1=0;s2=0;
fork=0:
（n-1）
x=a+h*k;
s1=s1+feval（fun,x）;
end
fork=0:
（n-1）
x=a+h*（k+1/2）;
s2=s2+feval（fun,x）;
end
s=h/6*（feval（fun,a）+feval（fun,b）+2*s1+4*s2）;
高斯求积公式求解：
functiong=gaosi（fname,a,b,n,m）
switchm
case1
t=0;
A=1;
case2
t=[-1/sqrt（3）,1/sqrt（3）];
A=[1,1];
case3
t=[-sqrt（0.6）,0.0,sqrt（0.6）];
A=[5/9,8/9,5/9];
case4
t=[-0.8612136,-0.339981,0.339981,0.861136];
A=[0.347855,0.652145,0.652145,0.347855];
case5
t=[-0.906180,-0.538469,0.0,0.538469,0.906180];
A=[0.236927,0.478629,0.568889,0.478629,0.236927];
case6
t=[-0.932470,-0.661209,-0.238619,0.238619,0.661209,0.932470];
A=[0.171325,0.360762,0.467914,0.467914,0.360762,0.171325];
otherwise
error
end
x=linspace（a,b,n+1）;
g=0;
fori=1:
n
g=g+gsint（fname,x（i）,x（i+1）,A,t）;
end
functiong=gsint（fname,a,b,A,t）
g=（b-a）/2*sum（A.*feval（fname,（b-a）/2*t+（a+b）/2））;
龙贝格求积公式求解：
function[t]=romberg（f,a,b,e）
t=zeros（15,4）;
t（1,1）=（b-a）/2*（f（a）+f（b））;
fork=2:
4
sum=0;
fori=1:
2^（k-2）
sum=sum+f（a+（2*i-1）*（b-a）/2^（k-1））;
end
t（k,1）=0.5*t（k-1,1）+（b-a）/2^（k-1）*sum;
fori=2:
k
t（k,i）=（4^（i-1）*t（k,i-1）-t（k-1,i-1））/（4^（i-1）-1）;
end
end
fork=5:
15
sum=0;
fori=1:
2^（k-2）
sum=sum+f（a+（2*i-1）*（b-a）/2^（k-1））;
end
t（k,1）=0.5*t（k-1,1）+（b-a）/2^（k-1）*sum;
fori=2:
4
t（k,i）=（4^（i-1）*t（k,i-1）-t（k-1,i-1））/（4^（i-1）-1）;
end
ifk>6
ifabs（t（k,4）-t（k-1,4））disp（[num2str（t（k,4））]）;
break;
end
end
end
ifk>=15
disp（['溢出']）;
end
3.3结果分析
复合辛普森公式
在命令窗口输入：
fun=inline（‘1./（1+x.^2）’）;
xinpusen（fun,01,10）
运行结果如下：
高斯求解高斯
在命令窗口输入：
gaosi（inline（‘1./（1+x.^2）’）,0,1,2,3）
运行结果如下：
龙贝格求积公式
运行结果如下：
3.4设计总结
实验结果显示每一个算法都接近真实值，但龙贝格算法相比较复合辛普森算法，高斯算法来说更加接近.对于代数精度来说，复合辛普森的代数精度为11，龙贝格代数精度为11，高斯代数精度为11。
可见代数精度相同时，龙贝格的求积精度最小，所以相同条件下龙贝格求积公式最能接近准确值。
这次课程设计，用MATLAB实验对数值积分进行了实现，用了3种不同的数值积分的方法，实现过程中发现了各种方法之间的区别和联系.并且在实验过程中，使自己对数值积分和MATLAB更加的熟悉.做到了学习和实践相联系。
设计题十一
非线性方程的数值解法的Matlab实现及其应用
4.1问题思路与设计思路
非线性是实际问题中经常出现的，并且在科学与工程计算中的地位越来越重要，很多我们熟悉的线性模型都是在一定条件下由非线性问题简化得到的。
非线性问题一般不存在直接的求解公式，故没有直接方法求解，都要使用迭代法求解，迭代法要求先给出根
的一个近似，若
且
，根据连续函数性质可知
在
内至少有一个实根，这时称
为方程
的有根区间，通常可通过逐次搜索法求得方程
的有根区间。
对于方程
，如果是线性函数，则它的求根是容易的。
牛顿法实质上是一种线性化方法，其基本思想是将非线性方程
逐步归结为某种线性方程来求解。
设已知方程
有近似根
，将函数
在点
展开，有
，于是方程
可近似地表示为：
，这是个线性方程
记其根为
，则
的计算公式为：
，
这就是牛顿法。
二分法的基本思想是将方程根的区间平分为两个小区间，把有根的小区间再平分为两个更小的区间，进一步考察根在哪个更小的区间内。
如此继续下去，直到求出满足精度要求的近似值。
4.2程序清单
迭代法：
求方程
的一个近似解，给定初值为
，误差不超过
。
function[p,k,err,P]=fixpt（f1021,p0,tol,max1）
P
（1）=p0;
fork=2:
max1
P（k）=feval（'f1021',P（k-1））;
k,err=abs（P（k）-P（k-1））
p=P（k）;
if（errbreak;
end
ifk==max1
disp（'maximumnumberofiterationsexceeded'）;
end
end
P=P;
定义一个名为f1021.m的函数文件
functiony=f1021（x）
y=sin（x）/x;
割线法：
求方程
的一个近似解，给定初值-1.5，-1.52，误差不超过
。
function[p1,err,k,y]=secant（f1042,p0,p1,delta,max1）
p0,p1,feval（'f1042',p0）,feval（'f1042',p1）,k=0;
fork=1:
max1
p2=p1-feval（'f1042',p1）*（p1-p0）/（feval（'f1042',p1）-feval（'f1042',p0））;
err=abs（p2-p1）;
p0=p1;
p1=p2;
p1,err,k,y=feval（'f1042',p1）;
if（errbreak,
end
end
先定义一个名为f1042.m的函数文件
functiony=f1042（x）
y=x^3-x+2;
建立一个主程序prog1042.m
clc
clear
secant（‘f1042’,-1.5,-1.52,10^（-6）,11）
牛顿法：
求解方程
的一个近似解，给定初值为1.2，误差不超过
。
function[p1,err,k,y]=newton（f1041,df1041,p0,delta,max1）
p0,feval（'f1041',p0）
fork=1:
max1
p1=p0-feval（'f1041',p0）/feval（'df1041',p0）;
err=abs（p1-p0）;
p0=p1;
p1,err,k,y=feval（'f1041',p1）
if（errbreak,
end
p1,err,k,y=feval（'f1041',p1）
end
先用m文件定义两个名为f1014.m和df1041.m的函数文件
functiony=f1041（x）
y=x^4-4*x+2;
functiony=df1041（x）
y=4*x^3-4;
建立一个主程序prog1041.m
clear
newton（‘f1041’,’df1041’,1.2,10^（-6）,18）
二分法：
functioner_fen（f,a,b,esp）;
f1=subs（f,a）;
f2=subs（f,b）;
iff1*f2>0
disp（'该方程在【a,b】上无解'）;
elseiff1==0
root=a;
elseiff2==0
root=b;
else
a0=a;
b0=b;
A=[];
whileabs（（b0-a0）/2）>=esp
half=（a0+b0）/2;
fa=subs（f,a0）;
fb=subs（f,b0）;
fhalf=subs（f,half）;
iffhalf==0
root=half;
break;
elseiffa*fhalf<0
b0=half;
else
a0=half;
end
A=[A,half];
end
root=（b0+a0）/2;
end
root
A
4.3结果分析
迭代法：
在命令窗口输入：
clc
clearall
fixpt（‘f1021’,0.5,10^（-5）,20）
运行结果如下：
k=2
err=0.4589
k=3
err=0.1052
k=4
err=0.0292
k=5
err=0.0078
k=6
err=0.0021
k=7
err=5.7408e-004
k=8
err=1.5525e-004
k=9
err=4.1975e-005
k=10
err=1.1350e-005
k=11
err=3.0688e-006
ans=0.8767
割线法：
在命令窗口输入：
clc
clear
secant（'f1042',-1.
升级会员 