数值分析综述报告.doc

《数值分析综述报告.doc》由会员分享，可在线阅读，更多相关《数值分析综述报告.doc（33页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、,淮 阴 工 学 院数值分析考试基于Matlab的方法综合应用报告班 级： 金融1121 姓 名： 姚婷婷 学 号： 1124104129 成 绩： 数 理 学 院 2014年6月7日数值分析课程综述报告前言：数值分析也称计算方法，它与计算工具的发展密切相关。数值分析是一门为科学计算提供必需的理论基础和有效、实用方法的数学课程，它的任务是研究求解各类数学问题的数值方法和有关的理论。 正文：第一章 近似计算与误差分析1、产生误差的原因：模型误差；观测误差；截断误差；舍入误差。2、四则运算的误差：加减法运算乘法运算 除法运算：3、科学表示法、有效数字、近似值的精度任何一个实数都可以表示成如下的形式

2、： 其中：是正整数，是整数，如果是数的近似值 并且 则称该近似值具有位有效数字（significant digit）。此时，该近似值的相对误差为 另一方面，若已知 那么，即：至少有位有效数字。例：取其近似值如下：x*=3.14x* =3.14159x*=3.1415x*=3.141第二章 线性方程组在科学计算中，问题的本身就是求解线性方程组，许多问题的求解需要最后也归结为线性方程组的求解，所以线性方程组的求解是科学计算中最常见的问题。对于线性方程组的求解一般有两种方法：（1) 直接法：高斯消去法；（2) 间接法：各种迭代法。(1) 高斯消去法：求解思路：先消元，即按一定的规律逐步消去未知量，将

3、方程组化为等价的上（或下）三角形方程组；然后进行回代，即由上三角形方程组逐个求出；高斯（列、全）主元素消去法,及在消元的每一步选取（列）主元素列中绝对值最大的元取做主元素，计算步骤：消元过程：按列选主元、行交换、消元计算；回代过程；高斯列主元素消去法的MATLAB 实现：。第三章 解线性方程组的迭代法通常逆矩阵不易求得，特别是对于大型的线性方程组，需要用迭代法求解。用迭代法求解线性方程组，要把线性方程组写成等价的形式，右边写为迭代格式，如：2、关于迭代法收敛性的两个重要结论：充分必要条件是：矩阵的谱半径；充分条件是：矩阵的某个算子范数。3、线性方程组的迭代法主要有Jacobian迭代法，Gau

4、ss-Seidel迭代法。Jacobian 迭代法: Gauss-Seidel 迭代法: (3.7)Jacobian 迭代法与G-S迭代法比较： （3.8）式（3.7）和 (3.8) 表明：Gauss-Seidel 迭代法在计算第次迭代的第个分量时，及时地利用了在此步迭代中得到的新的迭代值：，由于第步的迭代值通常比第步的迭代值更接近方程组的精确解，所以，在Jacobian迭代法和GS迭代法都收敛的情况下，Gauss-Seidel迭代法的收敛速度比Jacobian迭代法的收敛速度高。例题：用MATLAB函数normrdn生成5阶矩阵和向量分别构造线性方程组的Jacobi迭代格式和G-S迭代格式，

5、并判断收敛性。Jacobian迭代法和GS迭代法程序如下：clc;clear all;%1、生成M和bM=normrnd(1,2,5)b=normrnd(1,2,5,1)%Jacobian迭代法M1=D(L+U)f1=Dbrho=max(abs(eig(M1);R=1e-08; %设定的一个收敛标准switch sign(1-rho) case -1 disp(the Jocobian method is not applicable) otherwise x(:,1)=normrnd(0,9,5,1); k=1 while k=R k=k+1; else X=x(:,k+1); disp(J

6、acobian迭代法迭代次数为：) IterN=k %Jacobian迭代法迭代次数 break end endend%Causs-Seidel迭代法 M2=(D-L)Uf2=(D-L)brho=max(abs(eig(M2);R=1e-08;switch sign(1-rho) case -1 disp(the auss-seidel method is not applicable) otherwise x(:,1)=normrnd(0,9,5,1); k=1 while k=R k=k+1; else X=x(:,k+1) disp(Causs-Seidel迭代法迭代次数为：) Iter

7、N=k break end endend第四章 非线性代数方程（组）的数值解法：一、二分法：首先要确定适当的包含根的区间，这可以依据闭区间上连续函数的介值定理来确定，例如该方程：对于该方程所以该方程至少有一个实根位于区间，图像表明该区间中只含有一个实根；用表示方程在区间上的精确解，对于给定的精度要求，取区间的中点，并按下式进行判断： （1）以为例， 如果,没有达到精度要求，令，并重复（1）的迭代过程；如果，那么，必有，因为。即区间内的任何一点都可以作为方程的近根，特别地，可取做为近似解。二、牛顿迭代法：任取初始值上过点（）的切线方程为：与轴交于点，过点的切线方程为与轴交于点， ， ，如此下去得

8、牛顿迭代公式：例题：考虑如下三阶非线性方程组： 其中 取适当的迭代初值，用Newton迭代法求方程组的数值解程序：%Newton迭代法求解x=sym(x,clear);y=sym(y,clear);syms z; F=x2+y2+a2*z2/2-a2; x+y-a; (2*x-a)2+(2*y-a)2+a2*(z-b)/4; Fx=diff(F,x,1);Fy=diff(F,y,1);Fz=diff(F,z,1); DF=Fx Fy Fz; F=(x,y,z)x2+y2+a2*z2/2-a2;x+y-a; (2*x-a)2+(2*y-a)2+a2*(z-b)/4; %Newton迭代法求解过程

9、Fr=10-10; Er=10-10; x0=-a/4;-a;-a/2; x0=-a/4;-a;a/2;k=1; X(:,1)=x0;while norm(F(X(1,k),X(2,k),X(3,k)-0;0;0,2)=Fr tic; f(:,k)=F(X(1,k),X(2,k),X(3,k); J=subs(DF,x,X(1,k); J=subs(J,y,X(2,k);J=subs(J,z,X(3,k); X(:,k+1)=X(:,k)-Jf(:,k); t(k)=toc; if norm(X(:,k+1)-X(:,k),2)Er break end k=k+1;enddisp(Newton

10、迭代法结果为：);disp(X(:,end);运行结果：Newton迭代法结果为： 3.4291 0.46580.6535 第五章 插值一、插值：插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。二、常用差值法：拉格朗日（Lagrange）插值法、牛顿（Newton）插值法1、拉格朗日（Lagrange）插值法：拉格朗日插值多项式简单的证明因为拉格朗日插值多项式的基函数有如下的性质： 所以拉格朗日插值多项式 满足插值的条件。插值多项式拉格朗日插值法的不足在实际问题中，观测的数据可能会不断增加，如果用拉格朗日插值公式构造插值多项式，那么，每当增加

11、数据就要重新计算多项式的系数，由此增加许多不必要的计算工作量。2、（三次）样条（Spline）插值插值条件要求 要求所求的插值多项式（三次样条函数） 在每个区间，是次数不超过三次的多项式；，； 在区间.上具有二阶连续导数。例题：在夏季，较大湖泊的水体按深度被跃变层分为上部的变温层和下部的均温层。水体的分层化对环境工程中污染问题的研究具有重要的意义，例如，有机物的分解将导致被跃变层隔离的底部水体中氧急剧减少。按温度随水深的变化曲线，跃变层位于水深处： 现有美国普拉特湖（Platte Lake） 的一组数据：深度(): 02.34.99.113.718.322.927.2温度(): 22.822.

12、822.820.613.911.711.111.11. 试利用Lagrange差值方法求温度随水深（近似）变化函数表达式；2. 试利用三次样条差值方法（应用Matlab函数csape）求温度随水深（近似）变化函数表达式；3. 画出插值函数和曲线，并与原始插值数据图像作比较程序代码：函数文件程序：function Ln = my_Fun(x, XI, YI)if isa(x, sym) = 1; n = length(XI) - 1; Ln = 0; Pn = sym(ones(n + 1, 1); for k = 1 : n + 1 for i = 1 : n + 1 if i = k Pn(

13、k) = Pn(k)*(x - XI(i)/(XI(k) - XI(i); else Pn(k) = Pn(k); end end Ln = Ln + Pn(k)*YI(k); endelse n = length(XI) - 1; L = ones(n + 1, length(x); Ln = zeros(size(x); for k = 1 : n + 1 for i = 1 : n + 1 if i = k L(k,:) = L(k,:).*(x - XI(i)./(XI(k) - XI(i); else L(k,:) = L(k,:); end end Ln = Ln + L(k, :

14、).*YI(k); endEnd主文件程序：clc;clear all;close all;%问题1.用Lagrange差值方法求温度随水深（近似）变化图像t = sym(t, real);x=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2;T=22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1;X=linspace(0,27.2,275);Ln=my_Fun(X,x,T);figure(1)plot(x,T,r*,markersize,15) xlabel(深度x);ylabel(温度T);title(原始的散点图);pause(3)hold

15、on plot(X,Ln,b-,linewidth,3);xlabel(深度x);ylabel(温度T);title(Lagrange差值图);%求出温度随水深（近似）变化函数表达式I_poly = my_Fun(t, x, T);I_poly = simple(I_poly)I_poly = sym2poly(I_poly) %2.用三次样条差值方法求温度随水深（近似）变化函数 表达式figure(2)hold on y=interp1(x,T,X,spline);plot(X,y,b-,linewidth,3); %用三次样条插值方法画出图像xlabel(深度x);ylabel(温度T);

16、title(三次样条法差值图); PP=csape(x,T,2,2,0,0); Coefs = PP.coefs %3.两种差值函数图象比较figure(3)plot(X,Ln,-,color,1, 0, 0,LineWidth, 3)xlabel(深度x);ylabel(温度T);title(Lagrange差值图);hold on pause(3)fnplt(PP,b,3,0,28) %函数作图xlabel(深度x);ylabel(温度T);title(三次样条法差值图);pause(3)plot(x,T,-,linewidth,3,markersize,10)运行结果:I_poly =(

17、11156832321404503197312500000*t7 - 942692507666801544174224450000*t6 + 29858157926737187739430003070000*t5 - 433725749060135044183186809635000*t4 + 2850751947133625442244870404829250*t3 - 8285950923799759686915051372346805*t2 + 8477820147617239232740036583612328*t + 512657632901869980094703083834211

18、328)/22484983899204823688364170343605760I_poly = Columns 1 through 6 0.0000 -0.0000 0.0013 -0.0193 0.1268 -0.3685 Columns 7 through 8 0.3770 22.8000Coefs = 0.0022 -0.0000 -0.0115 22.8000 -0.0092 0.0150 0.0230 22.8000 -0.0114 -0.0565 -0.0850 22.8000 0.0297 -0.2004 -1.1640 20.6000 -0.0153 0.2099 -1.12

19、00 13.9000 0.0017 -0.0014 -0.1605 11.7000 -0.0017 0.0223 -0.0640 11.1000 图示如下： 第六章 最小二乘拟合与最佳逼近一、最小二乘拟合加权最小二乘法逼近准则：最小二乘逼近多项式必须满足如下必要条件： 即满足法方程组: 例题：下面是一处地质岩层断面高程（深度）的测量数据。水平距离():00.201.002.103.505.006.807.509.0011.212.0高度():h 1.641.581.681.841.580.860.390.310.390.770.86试利用最小二乘法求满足 即误差不超过的最低次数的拟合多项式,写

20、出该多项式的表达式;画出数据散点图和该多项式曲线.程序：clc;clear all;close all;x=linspace(0,12,12);x=0 200 1000 2100 3500 5000 6800 7500 9000 11200 12000;h=1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86;plot(x,h,*,markersize,8) %画出给出数据的散点图xlabel(水平距离x);ylabel(高度h);title(岩层断面水平距离x和高度h的散点图);figure(1) %1.最小二乘法拟合 n=3;P, S

21、 = polyfit(x, h, n)Pm = polyval(P, x)R = sqrt(sum(h - Pm).2) %误差 t = linspace(x(1), x(end), 12);Poly = polyval(P, t);figure(2)plot(x, h, ro, t, Poly, LineWidth, 3, markersize, 8)set(gca,FontSize,18)legend(The Data, The Fitting Curve, 1)title(Curve Fitting by Least Square Approximation, fontsize, 18)

22、第七章 微积分的数值方法一、数值微分如果给定函数 的关系式比较复杂或者 未知，而仅仅知道在个相异点，处的函数值，那么，我们可以利用函数的插值多项式的导数作为函数导数的近似值 因而有 这里需要说明一点的是，尽管和的函数值可能相差不多，但是它们的导数有可能相差很大。二、数值积分数值积分方法最大的优点是将复杂的函数积分转化为相对简单的加、减、乘、除运算。对于给定的函数，如果 的函数关系式比较复杂或.未知，而仅仅知道该函数在个点处的函数值，那么可以利用函数的插值多项式的积分作为函数在上的定积分的近似值。1、Newton-Cotes求积公式根据Lagrange插值多项式有 令 得Newton-Cotes

23、求积公式特别地，当取插值节点为时有两点公式（梯形公式）： 三点公式（Simpson公式）： 例题：下面是一处地质岩层断面上部边缘的深度测量数据。水平距离():00.201.002.103.505.006.807.509.0011.212.0深度(): 1.641.581.681.841.580.860.390.310.390.770.86表1试利用复化的梯形求积法求该组数据所在曲线与基准线（轴）在范围内所围成图形面积.画出数据散点图和图形的示意图.试利用复化的Simpson求积法求该组数据所在曲线与基准线（轴）在范围内所围成图形面积.画出数据散点图和图形的示意图.解答如下：程序：% 题目给出的

24、数据XI = 0 0.20 1.00 2.10 3.50 5.00 6.80 7.50 9.00 11.2 12.00;YI = -1.64 1.58 1.68 1.84 1.58 0.86 0.39 0.31 0.39 0.77 0.86;n = length(XI) - 1; %此时n=10M = 60;x = linspace(XI(1), XI(end), M +1); %设置线性空间% 问题1.1 T1 = 0; %初值%将图像分成若干个小区间，在每个区间内，求小梯形的面积，并累加起来，就是题目所要求的面积for k = 2 : n + 1 T1 = T1 + (YI(k) + YI

25、(k-1)*(XI(k) - XI(k-1)/2;enddisp(复化的梯形求积法求得的图形面积：)T_Trape = T1 %累加后的图形面积figure(1) %画出函数图像set(gca,FontSize, 20)patch(XI(1), XI, XI(end), 0, YI, 0, 0.5 1 0.5) %颜色的填充hold onplot(XI, YI, o-, XI, 0*XI, k, LineWidth, 3, markersize, 8)title(复化的梯形求积法求围成图形面积)xlabel(水平距离x)ylabel(深度y)运行结果：复化的梯形求积法求得的图形面积：T_Tra

26、pe = -11.6090图像：程序：函数M文件：function Ln = Lagrange_Fun_01(x, XI, YI)if isa(x, sym) = 1; n = length(XI) - 1; Ln = 0; Pn = sym(ones(n + 1, 1); for k = 1 : n + 1 for i = 1 : n + 1 if i = k Pn(k) = Pn(k)*(x - XI(i)/(XI(k) - XI(i); else Pn(k) = Pn(k); end end Ln = Ln + Pn(k)*YI(k); endelse n = length(XI) -

27、1; L = ones(n + 1, length(x); Ln = zeros(size(x); for k = 1 : n + 1 for i = 1 : n + 1 if i = k L(k,:) = L(k,:).*(x - XI(i)./(XI(k) - XI(i); else L(k,:) = L(k,:); end end Ln = Ln + L(k, :).*YI(k); endend主程序：% 问题1.2 和 1.3%问题3是问题2的延伸S2 = 0;%Simpson复化求积法，应用同样的面积累加，得到最终的图形面积for k = 1: 2: n - 1 S2 = S2 +

28、(XI(k+2) - XI(k)*(YI(k) + 4*YI(k+1) + YI(k+2)/6; X(k, :) = linspace(XI(k), XI(k+2), 20); Y(k, :) = Lagrange_Fun_01(X(k,:), XI(k : k+2), YI(k : k+2);enddisp(Simpson复化求积法求得的图形面积：)S_Simpson = S2 %Simpson复化求积法求得的面积figure(2)set(gca,FontSize, 20)hold onfor k = 1: 2 : n-1 patch(XI(k), X(k,:), XI(k +2), 0,

29、Y(k, :), 0, 0.5 1 0.5) plot(X(k, :), Y(k, :), LineWidth, 3)endplot(XI, YI, o, XI, 0*XI, k, LineWidth, 3, markersize, 8)title(复化的Simpson求积法求围成图形面积)xlabel(水平距离x)ylabel(深度y)运行结果：Simpson复化求积法求得的图形面积：S_Simpson = -11.9128图像：第八章 微分方程（组）初值问题数值方法引言：微分方程（组）是科学研究和工程应用中最常用的数学模型之一。但是基绝大多数的常微分方程和常微分方程组得不到解析解和实际应用

30、中往往只需要知道常微分方程（组）的解在某些点处的函数值这两点原因，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。一、 Euler 方法最简单的数值解法1、Euler 法假设要求在点,， 处初值问题（4）的解的近似值。首先对式（4）的两端积分，得 对于该式的右边，如果被积函数用积分下限处的函数值代替被积函数作积分，则有进而得到下式给出的递推算法Euler 方法2、Euler 方法的改进改进的Euler 方法Euler 方法的精度高，其原因在于：在推导Euler 方法时，我们是用待求解函数在一点处的

31、变化率 （1）代替在区间上的平均变化率；而在推导改进的Euler 方法时，我们是用待求解函数在两点处变化率的平均值（2）代替在区间上的平均变化率；显然，通常式（2）比式（1）更接近于在区间上的平均变化率。由此，人们受到启发：适当地选取区间上函数多个点处的变化率，用它们加权平均值代替在区间上的平均变化率，近似解的精度应更高。二、 RungeKutta法：二阶的RungeKutta法三阶的RungeKutta法四阶的RungeKutta法RungeKutta法是按选取区间上函数变化率数目的多少和截断误差的阶数来区分的一系列方法。例题：给定如下常微分方程Cauchy 问题： 分别取参数，初值和步长：

32、 ，。 取，在区间上分别利用Euler 法，4阶R-K公式和matlab的ode45求问题（1）的数值解；程序：clcclear alla=1;b=1;F = (t, y) -a*y-t2*exp(-0.5*t)*sin(t)+b;%欧拉法、matlab的ode45方法求解及画图y0 = -1;Tf = 2*pi;Tspan = 0, Tf;Options = odeset(RelTol, 1e-5, AbsTol, 1e-6);T, z = ode45(F, Tspan, y0,Options);h = 0.2;t1 = 0:h:Tf;Lt = length(t1);Y(1) = y0;fo

33、r k = 2 : Lt Y(k) = Y(k-1) + h*F(t1(k-1), Y(k-1); endset(gca,Fontsize, 24)plot(T, z, t1, Y, t1, Y, o, LineWidth, 4, markersize, 5) %4阶R-K公式求解及画图tt = 0 : h : Tf;c = 1;YY(1) = y0;for k = 1: length(tt) -1; K1(k) = F(tt(k), YY(k); z(k) = YY(k) + 2h*K1(k); K2(k) = F(tt(k)+h/2, z(k); z(k) = YY(k) + 2h*K2(

34、k); K3(k) = F(tt(k)+h/2, z(k); z(k) = YY(:,k) + h*K3(k); K4(k) = F(tt(k)+h, z(k); YY(k+1) = YY(k) + 6h*(K1(k) + 2*K2(k) + 2*K3(k) + K4(k);endhold onplot(tt, YY, LineWidth, 3)title(An Example of Applying The Fourth-Order Rong - Kutta Method)legend(Euler Curve, The Curve Obtained By Order 4 R-K,The Curve Obtained By ode45)xlabel(Time)ylabel(Orbit)%输出三种方法最后的数值解C=cell(2,3);C=Euler法,4阶R-K法,ode45法;Y(end),YY(end),z(end)运行结果：最终得到的三种方法的数值解分别为：C = Euler法 4阶R-K法 ode45法 2.1863 2.0838 2.0029参考文献：1、盖如栋，Lecture 1_引言、近似计算与误差分析，电子讲义，淮阴工学院，2014；2、盖如栋，数值分析实验，电子讲义，淮阴工学院，2013。