数值分析课程设计含代码

1、 成 绩 评 定 表学生姓名班级学号专 业信息与计算科学课程设计题目数值分析算法案例评语组长签字：成绩日期 20 年 月 日课程设计任务书学 院理学院专 业信息与计算科学学生姓名班级学号课程设计题目数值分析算法案例实践教学要求与任务:要求：格式以学校毕业论文格式要求为准，不准粘贴图片，尤其公式。对每个试验，要求有：实验基本原理，实验目的，实验内容及数据来源和实验结论。以班级为单位统一装订封皮。6月25日，十八周的周二交论文每人至少四个实验，最少15页任务（实验项目）：线性方程组数值解法 参考题目：（1) 列主元Gauss消去法（2）LU分解法插值法和数据拟合 参考题目：（1）Lagrange插

2、值（2）Newton插值（3）最小二乘拟合数值积分 参考题目：（1）复化Simpon积分（2）变步长的梯形积分公式（3）龙贝格求积公式常微分方程数值解 Runge-Kutta方法数值方法实际应用 用数值方法解决实际问题（自选）工作计划与进度安排:线性方程组数值解法 （4学时）插值法和数据拟合 （4学时）数值积分 （4学时）常微分方程数值解 （4学时）数值方法实际应用 （4学时）答辩 （4学时）指导教师： 201 年 月 日专业负责人：201 年 月 日学院教学副院长：201 年 月 日摘 要实验方法与理论方法是推动科学技术发展的两大基本方法，但有局限性。许多研究对象，由于空间或时间的限制，既不

3、可能用理论精确描述，也不能用实验手段实现。 数值模拟或称为科学计算突破了实验和理论科学的局限，在科技发展中起到越来越重要的作用。可以认为，科学计算已于实验、理论一起成为科学方法上不可或缺的三个主要手段。 计算数学的研究是科学计算的主要组成部分，而数值分析则是计算数学的核心。数值计算是研究使用计算机来解决各种数学问题的近似计算方法与理论，其任务是提供在计算机上可解的、理论可靠的、计算复杂性低的各种常用算法。数值分析的主要内容：1）、数值代数：求解线性和非线性方程组的解，分直接方法和间接方法两大类；2）、插值、曲线拟合和数值逼近；3）、数值微分和数值积分；4）、常微分和偏微分方程数值解法。 本文主

4、要通过Matlab软件，对数值分析中的一些问题进行求解，如列主元Gauss消去法，Lagrange插值多项式，复化Simpson公式，Runge-Kutta方法以及数值分析在实际问题中的应用，并在求解的过程中更加熟识这门课程的主要内容，以及加强对课程知识的掌握。在学习与设计计算方法时，从数学理论角度，学会分析方法的误差、收敛性和稳定性，保证计算方法的准确性；从实际应用的角度出发，掌握计算方法的结构与流程，能够把计算方法转换为可在计算机上直接处理的程序，保证算法的可用性。关键词：列主元Gauss消去法；Lagrange插值；复化Simpson公式；Runge-Kutta目 录实验一 列主元Gau

5、ss消去法11.1 实验目的11.2 基本原理11.3 实验内容21.4 实验结论3实验二 拉格朗日插值多项式42.1 实验目的42.2 基本原理42.3 实验内容42.4 实验结论9实验三 复化Simpson求积公式103.1 实验目的103.2 基本原理103.3 实验内容103.4 实验结论12实验四 龙格-库塔(Runge-Kutta)方法134.1 实验目的134.2 基本原理134.3 实验内容144.4 实验结论15实验五 数值方法实际应用165.1 实验目的165.2 基本原理165.3 实验内容165.4 实验结论22参考文献23实验一 列主元Gauss消去法1.1 实验目的

6、1) 理解列主元消去法的原理；2) 熟悉列主元消去法的计算步骤，能用代码编写；3) 解决实际问题。1.2 基本原理在顺序Gauss消去法中，必须要求；否则无法进行计算。即使，但其绝对值很小，由于舍入误差的影响，也可能会引起很大的误差，从而使上述方法失效。为了使消元过程中减小舍入误差和不至于中断，可以按照不同的自然顺序进行消元。在第k步消元时，增广矩阵为 （1.1）不一定选取作为主元，而从同列中选取绝对值最大的作为主元素，即 （1.2） 若，此时矩阵不可逆，方程的解不确定，则停止计算；否则，当r>k时，则其增广矩阵中交换第k行和第r行，即 （1.3）使成为主元。然后再按Gauss消去法进行

7、消元运算。于是就得到列主元Gauss消去法。1.3 实验内容1.3.1 程序来源 首先建立一个gaussMethod.m的文件，用来实现列主元的消去方法。文件内容如下：function x=gaussMethod(A,b)%高斯列主元消去法，要求系数矩阵非奇异的n = size(A,1);if abs(det(A)<= 1e-8 error('系数矩阵是奇异的'); return;endfor k=1:n ak = max(abs(A(k:n,k); index = find(A(:,k)=ak); if length(index) = 0 index = find(A(

8、:,k)=-ak); end %交换列主元 temp = A(index,:); A(index,:) = A(k,:); A(k,:) = temp; temp = b(index);b(index) = b(k); b(k) = temp; %消元过程 for i=k+1:n m=A(i,k)/A(k,k); %消除列元素 A(i,k+1:n)=A(i,k+1:n)-m*A(k,k+1:n); b(i)=b(i)-m*b(k); endend %回代过程 x(n)=b(n)/A(n,n); for k=n-1:-1:1; x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)'

9、;)/A(k,k); end; end然后调用gaussMethod函数，来实现列主元的高斯消去法。建立一个文件gauss，内容如下：cleardisp('*')x=gaussMethod(input('请输入系数矩阵：'),input('请输入常数列：')disp('*')1.3.2 实例分析例：在Matlab上，利用列主元法求线性方程组的解: 解：运行程序，按照提示输入方程的系数矩阵及常数列，如下所示：*请输入系数矩阵：1 2 1 4;2 0 4 3;4 2 2 1;-3 1 3 2请输入常数列：13;28;20;6x = 3

10、 -1 4 2* 即该方程的解为：1.4 实验结论把向量计算得到的解带入方程组，验证正确性，和其他的方法比较，列主元具有一定的简单性，比较容易实现。避免使用其他方法的误差或不能进行性。而列主元也有一定的限制，要求行列式的值不为0。实验二 拉格朗日插值多项式2.1 实验目的1）熟悉简单的拉格朗日插值多项式的基本概念；2）熟悉Lagrange公式及源代码，会利用它来计算基本函数；3）能构造出正确的插值多项式；2.2 基本原理设函数在区间a,b上有定义，且已知在点上的值若存在一个次数不超过n的多项式 （2.1）使其满足 （2.2）则称为的n次插值多项式，称点为插值节点，称条件（2.2）为插值条件。包

11、含插值节点的区间成为插值区间。通过平面上不同的两点可以确定一条直线经过这两点，就是拉格朗日线性插值问题，对于不在同一直线的三点得到的插值多项式为抛物线。拉格朗日是比较基础的方法，本身比较容易实现，容易理解。给定n+1个不同节点，构造的n次拉格朗日插值多项式： ， （2.3）2.3 实验内容2.3.1 程序来源首先建立一个Lagrange.m的文件，用来实现Lagrange插值。文件内容如下：%输入：x是插值节点横坐标向量；y是插值节点对应纵坐标向量%输出：C是拉格朗日插值多项式的系数矩阵；L是插值函数系数矩阵functionC,L=Lagrange(x,y)w=length(x);n=w-1;

12、L=zeros(w,w);for k=1:n+1 V=1; for j=1:n+1 if k=j V=conv(V,poly(x(j)/(x(k)-x(j); end end L(k,:)=V;endC=y*L然后调用Lagrange函数，来实现Lagrange插值法。建立一个文件Lg，内容如下：cleardisp('*')x=input('请输入已知点的横坐标组：');y=input('请输入已知点的纵坐标：');C,L=Lagrange(x,y);yi=polyval(C,input('请输入需要计算得横坐标组：')xx=1.

13、5:0.05:6.5;yy=polyval(C,xx);plot(xx,yy,x,y,'o')disp('*')2.3.2 实例分析例1 有4对数据（1.6,3.3），（2.7,4.22），（3.9,5.61），（5.6,2.94），写出这4个数据点的Lagrange插值公式，并计算出横坐标组xi=2.101,4.234时对应的纵坐标值。解：4个数据点的Lagrange插值公式为： 运行程序，按照提示输入已知点的横坐标组、纵坐标组及需要计算得横坐标组，如下所示：*请输入已知点的横坐标组：1.6,2.7,3.9,5.6请输入已知点的纵坐标：3.3,1.22,5.6

14、1,2.94C = -1.0539 11.0551 -34.4933 34.5053请输入需要计算得横坐标组：2.101,4.234yi = 1.0596 6.6457*即 输出图形：图2.1 输出拟合曲线例2 将区间-5，5等分5份、10份，求函数拉格朗日差值多项式，做出函数原图像，观察龙格现象。解：首先将区间五等分，取各端点坐标拟合曲线，输入：clearx=-5:2:5;y=1./(1+x.2);C,L=Lagrange(x,y);输出结果：C = 0.0000 0.0019 -0.0000 -0.0692 0.0000 0.5673即 然后输出图形，比较拟合效果。在matlab中输入：x

15、x=-5:0.1:5;yy=polyval(C,xx);hold onplot(xx,yy,'-',x,y,'o')xp=-5:0.01:5;z=1./(1+xp.2);plot(xp,z,'r')legend('拉格朗日插值曲线','插值点','原曲线')输出五等分插值龙格现象图形：图2.2 五等分插值龙格现象图形 十等分的过程与五等分基本相似，如下输入：clear,clfx=-5:1:5;y=1./(1+x.2);C,L=Lagrange(x,y);输出：C = -0.0000 -0.0000

16、0.0013 -0.0000 -0.0244 0.0000 0.1974 -0.0000 -0.6742 -0.0000 1.0000即 然后输出图形，比较拟合效果。在matlab中输入：xx=-5:0.1:5;yy=polyval(C,xx);hold onplot(xx,yy,'-',x,y,'o')xp=-5:0.01:5;z=1./(1+xp.2);plot(xp,z,'r') legend('拉格朗日插值曲线','插值点','原曲线')输出十等分插值图形：图2.3 十等分插值龙格现象图形2

17、.4 实验结论通过图像可以得出：1）并不是插值节点越多，插值多项式逼近函数效果就越好；2）误差较大的地方，是在插值区间两端点附近出现；3）求插值多项式，不需要求解线性方程组，当数据比较多时，此公式才能 显示出优越性；4）函数值可以用符号形式表示，数据点未确定的纵坐标可以用多项式表示；实验三 复化Simpson求积公式3.1 实验目的1）了解单独的Simpson公式原理及几何意义；2）学会复化Simpson求积公式的编程与应用；3）掌握Mat lab提供的计算积分的各种函数的使用方法；3.2 基本原理 复化Simpson公式是一种比较实用的积分方法，可以给出误差估计。首先将区间a,b N等分，子

18、区间的长度为 （3.1）在每个子区间上采用Simpson公式，在用Simpson公式时，还要将子区间再二等分，因此有2N+1个分点。即 （3.2）经推导得到， （3.3）称为为复化Simpson值，称（3.3）式为复化Simpson公式。3.3 实验内容3.3.1 程序来源 编写复化Simpson求积函数（函数名：s_quad.m）Function I=S_quad(x,y)% 复化求积公式% x为被积函数自变量的等距节点；y为被积函数在节点处的函数值。n=length(x);m=length(y); % 积分自变量的节点数应与它的函数值个数相同；if n=m error ('The

19、length of X and Y must be equal'); return;endif rem(n-1,2)=0 % 如果n-1不能被2整除，则调用复化公式 error ('节点数不满足要求'); return;endN=(n-1)/2;h=(x(n)-x(1)/N;a=zeros (1,n);for k=1:Na(2*k-1)=a(2*k-1)+1;a(2*k)=a(2*k)+4; a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);然后调用s_quad函数，来实现复化Simpson公式法。建立一个文件SPS，内容如下：cleardi

20、sp('*')x=input('请输入积分上下限及点间的间隔（例如-1:0.1:1）：');y=input('请输入被积公式：y=');I=S_quad(x,y);disp('得出积分值I=')disp(I);disp('*')3.3.2 实例分析：例1 用复化Simpson公式求积分，在积分区间中点与点之间的间隔取为 0.1。解：运行程序，按照提示输入积分上下限、点间的间隔及被积公式，如下所示：*请输入积分上下限及点间的间隔（例如-1:0.1:1）：-1:0.1:1请输入被积公式：y=exp(-x.2)得出积分值

21、I= 1.4936* 真值为：1.4937例2 计算积分，将区间8等分。解：运行程序，按照提示输入积分上下限、等分后的区间长度及被积公式，如下 所示：*请输入积分上下限及点间的间隔（例如-1:0.1:1）：0:0.125:1请输入被积公式：y=x./(4+x.2)得出积分值I= 0.1116* 真值为：0.1115723.4 实验结论复化Simpson计算所得的结果误差较小，精度较高，更适合科学计算与应用，且公式具有收敛性，稳定性良好。实验四 龙格-库塔(Runge-Kutta)方法4.1 实验目的1）熟悉Runge-Kutta常微分方程初值问题的基本原理2）了解Runge-Kutta常微分方

22、程初值问题的计算3）能编程实现Runge-Kutta常微分方程初值问题4.2 基本原理 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有： (4.1)当用点处的斜率近似值与右端点处的斜率的算术平均值作为平均斜率的近似值，那么就会得到二阶精度的改进欧拉公式： (4.2)依次类推，如果在区间内多预估几个点上的斜率值，并用他们的加权平均数作为平均斜率的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得四阶龙格库塔公式，也就

23、是在工程中应用广泛的经典龙格库塔算法。 四阶龙格库塔方法基本公式： (4.3)4.3 实验内容4.3.1 程序来源写出四阶Runge-Kutta方法源代码的得.m文件：function R=rk4(f,a,b,ya,M)% a,b分别为左右端点，ya为初始值，M为步长幅数。h=(b-a)/M;T=zeros(1,M+1);Y=zeros(1,M+1);T=a:h:b;Y(1)=ya;for j=1:M k1=feval(f,T(j),Y(j); k2=feval(f,T(j)+h/2,Y(j)+k1*h/2); k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);k4=feval

24、(f,T(j)+h,Y(j)+h*k3); Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endR=T' Y'4.3.2 实例分析：例 用四阶龙格库塔方法求解初值问题取步长h=0.1。并将a,b,ya,M分别取定0,1,0,10。解：首先写出执行函数的tt.m文件，其中tt表示f，即为求导函数：function f=tt(t,y)f=t2+t-y;再输入调用四阶Runge-Kutta方法，执行语句如下：>> R=rk4('tt',0,1,0,10)输出：R = 0 0 0.30000000000000 0.3075172207808

25、9 0.40000000000000 0.46666197888861 0.50000000000000 0.64581185656207 0.60000000000000 0.84496198189452 0.80000000000000 1.30326225710000 0.90000000000000 1.562412395020554.4 实验结论龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算在一些点上的值，如四阶龙格-库塔方法每计算一步需要计算四次的值，这给实际计算带来一定的复杂性。实验五 数值方法实际应用5.1 实验

26、目的1）了解梯形公式原理及几何意义；2）学会复化梯形求积公式的编程与应用；3）掌握Matlab提供的计算积分的各种函数的使用方法；5.2 基本原理 在实际问题中，往往会遇到一些困难。有些函数找不到用初等函数表示的原函数，例如，对于积分 (5.1)而言，不存在用初等函数表示的原函数。而有些函数虽然能找到原函数，但计算过于复杂，例如，椭圆型积分 (5.2) 而有些情况下，只能知道某些点处的函数值，并没有函数的具体表达式。这些情况，使我们有必要研究积分的数值计算问题。下面我们就以梯形公式为例做以说明。 所谓梯形求积公式就是用梯形面积来近似曲边梯形面积，利用梯形公式和连续增加a,b的区间数来逼近： (

27、5.3)第j次循环在个等距节点处对采样。5.3 实验内容 卫星轨道是一个椭圆，椭圆周长计算公式是，这里a是椭圆半长轴，c是地球中心与轨道中心（椭圆中心）的距离，记h为近地点距离，H为远地点距离，R=637km为地球半径，则我国第一颗人造卫星近地点距离h=439km，远地点距离H=2384，试求卫星轨道的周长。解: 第一步：先利用Matlab画出被积函数的图形。 输入程序如下：clearH=2384;h=439;R=6371;a=(2*R+H+h)/2c=(H-h)/2x=0:0.1:pi/2;y=sqrt(1-(c/a)2*(sin(x).2);plot(x,y,'-')tit

28、le('梯形法则');xlabel('x');ylabel('y');输出结果：a = 7.782500000000000e+003c = 9.725000000000000e+002输出图形：图5.1 被积函数的图形 第二步：应用数值积分梯形公式。 首先建立一个名为trapezg.m的M文件，程序如下：function I=trapezg(f_name3,a,b,n)format long%输出用15位数字表示n=n;h=(b-a)/n; x=a+(0:n)*h;f=feval(f_name3,x);I=h/2*(f(1)+f(n+1);if

29、n>1 I=I+h*sum(f(2:n);endh1=(b-a)/100;xc=a+(0:100)*h1; fc=feval(f_name3,xc);plot(xc,fc,'r');hold onxlabel('x');ylabel('y');plot(x,f)title('数值积分梯形效果图');plot(x,zeros(size(x),'.')for i=1:n;plot(x(i),x(i),0,f(i),end 然后建立一个名为f_name3.m的M文件定义函数，Matlab命令如下：function y=f_name3(x)y=sqrt(1-(9.725000e+002/