数值分析综述报告

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

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

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

4、(3.7)Jacobian 迭代法与G-S迭代法比拟：3.8式3.7和 (3.8) 说明：Gauss-Seidel 迭代法在计算第次迭代的第个分量时，及时地利用了在此步迭代中得到的新的迭代值：，由于第步的迭代值通常比第步的迭代值更接近方程组的准确解，所以，在Jacobian迭代法和GS迭代法都收敛的情况下，Gauss-Seidel迭代法的收敛速度比Jacobian迭代法的收敛速度高。例题：用MATLAB函数normrdn生成5阶矩阵和向量分别构造线性方程组的Jacobi迭代格式和G-S迭代格式，并判断收敛性。Jacobian迭代法和GS迭代法程序如下：clc;clear all;%1MbM=n

5、ormrnd(1,2,5)b=normrnd(1,2,5,1)%JacobianM1=D(L+U)f1=Dbrho=ma*(abs(eig(M1);R=1e-08; %switch sign(1-rho) case -1 disp(the Jocobian method is not applicable) otherwise *(:,1)=normrnd(0,9,5,1); k=1 while k=R k=k+1; else *=*(:,k+1); disp(Jacobian) IterN=k %Jacobian break end endend%Causs-Seidel M2=(D-L)U

6、f2=(D-L)brho=ma*(abs(eig(M2);R=1e-08;switch sign(1-rho) case -1 disp(the auss-seidel method is not applicable) otherwise *(:,1)=normrnd(0,9,5,1); k=1 while k=R k=k+1; else *=*(:,k+1) disp(Causs-Seidel) IterN=k break end endend第四章非线性代数方程组的数值解法：一、二分法：首先要确定适当的包含根的区间，这可以依据闭区间上连续函数的介值定理来确定，例如该方程：对于该方程所以该

7、方程至少有一个实根位于区间，图像说明该区间中只含有一个实根；用表示方程在区间上的准确解，对于给定的精度要求，取区间的中点，并按下式进展判断：1以为例，如果,没有到达精度要求，令，并重复1的迭代过程；如果，则，必有，因为。即区间内的任何一点都可以作为方程的近根，特别地，可取做为近似解。二、牛顿迭代法：任取初始值上过点的切线方程为：与轴交于点，过点的切线方程为与轴交于点，如此下去得牛顿迭代公式：例题：考虑如下三阶非线性方程组：其中取适当的迭代初值，用Newton迭代法求方程组的数值解程序：%Newton*=sym(*,clear);y=sym(y,clear);syms z; F=*2+y2+a2

8、*z2/2-a2; *+y-a; (2*-a)2+(2*y-a)2+a2*(z-b)/4; F*=diff(F,*,1);Fy=diff(F,y,1);Fz=diff(F,z,1); DF=F* Fy Fz; F=(*,y,z)*2+y2+a2*z2/2-a2;*+y-a; (2*-a)2+(2*y-a)2+a2*(z-b)/4; %NewtonFr=10-10; Er=10-10; *0=-a/4;-a;-a/2; *0=-a/4;-a;a/2;k=1; *(:,1)=*0;while norm(F(*(1,k),*(2,k),*(3,k)-0;0;0,2)=Fr tic; f(:,k)=F

9、(*(1,k),*(2,k),*(3,k); J=subs(DF,*,*(1,k); J=subs(J,y,*(2,k);J=subs(J,z,*(3,k); *(:,k+1)=*(:,k)-Jf(:,k); t(k)=toc; if norm(*(:,k+1)-*(:,k),2)Er breakend k=k+1;enddisp(Newton);disp(*(:,end);运行结果：Newton迭代法结果为：3.4291 0.46580.6535 第五章插值一、插值：插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。二、常用差值法：拉格朗日L

10、agrange插值法、牛顿Newton插值法1、拉格朗日Lagrange插值法：拉格朗日插值多项式 = 1 * GB2 简单的证明因为拉格朗日插值多项式的基函数有如下的性质：所以拉格朗日插值多项式满足插值的条件。插值多项式 = 3 * GB2 拉格朗日插值法的缺乏在实际问题中，观测的数据可能会不断增加，如果用拉格朗日插值公式构造插值多项式，则，每当增加数据就要重新计算多项式的系数，由此增加许多不必要的计算工作量。2、三次样条Spline插值 = 1 * GB2 插值条件要求要求所求的插值多项式三次样条函数 = 1 * GB3 在每个区间，是次数不超过三次的多项式； = 2 * GB3 ，； =

11、 3 * GB3 在区间.上具有二阶连续导数。例题：在夏季，较大湖泊的水体按深度被跃变层分为上部的变温层和下部的均温层。水体的分层化对环境工程中污染问题的研究具有重要的意义，例如，有机物的分解将导致被跃变层隔离的底部水体中氧急剧减少。按温度随水深的变化曲线，跃变层位于水深处：现有美国普拉特湖Platte Lake的一组数据：深度(): 02.34.99.113.718.322.927.2温度(): 22.822.822.820.613.911.711.111.11. 试利用Lagrange差值方法求温度随水深近似变化函数表达式；2. 试利用三次样条差值方法应用Matlab函数csape求温度随

12、水深近似变化函数表达式；3. 画出插值函数和曲线，并与原始插值数据图像作比拟程序代码：函数文件程序：function Ln = my_Fun(*, *I, YI)if isa(*, sym) = 1; n = length(*I) - 1; Ln = 0; Pn = sym(ones(n + 1, 1);for k = 1 : n + 1for i = 1 : n + 1if i = k Pn(k) = Pn(k)*(* - *I(i)/(*I(k) - *I(i);else Pn(k) = Pn(k);endend Ln = Ln + Pn(k)*YI(k);endelse n = leng

13、th(*I) - 1; L = ones(n + 1, length(*); Ln = zeros(size(*);for k = 1 : n + 1for i = 1 : n + 1if i = k L(k,:) = L(k,:).*(* - *I(i)./(*I(k) - *I(i);else L(k,:) = L(k,:);endend Ln = Ln + L(k, :).*YI(k);endEnd主文件程序：clc;clear all;close all;%1.Lagranget = sym(t, real);*=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2;T=

14、22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1;*=linspace(0,27.2,275);Ln=my_Fun(*,*,T);figure(1)plot(*,T,r*,markersize,15) *label(*);ylabel(T);title();pause(3)hold onplot(*,Ln,b-,linewidth,3);*label(*);ylabel(T);title(Lagrange);%I_poly = my_Fun(t, *, T);I_poly = simple(I_poly)I_poly = sym2poly(I_poly)%2. f

15、igure(2)hold ony=interp1(*,T,*,spline);plot(*,y,b-,linewidth,3); %*label(*);ylabel(T);title();PP=csape(*,T,2,2,0,0); Coefs = PP.coefs%3.figure(3)plot(*,Ln,-,color,1, 0, 0,LineWidth, 3)*label(*);ylabel(T);title(Lagrange);hold onpause(3)fnplt(PP,b,3,0,28) %*label(*);ylabel(T);title();pause(3)plot(*,T,

16、-,linewidth,3,markersize,10)运行结果:I_poly =(0*t7 - 00*t6 + 0000*t5 - 35000*t4 + *t3 - *t2 + *t + 28)/0I_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

17、 22.8000 0.0297 -0.2004 -1.1640 20.6000 -0.0153 0.2099 -1.1200 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

18、.860.390.310.390.770.86试利用最小二乘法求满足即误差不超过的最低次数的拟合多项式,写出该多项式的表达式;画出数据散点图和该多项式曲线.程序：clc;clear all;close all;*=linspace(0,12,12);*=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(*,h,*,markersize,8) %画出给出数据的散点图*label(*);ylabel(h);title(*

19、h);figure(1)%1. n=3;P, S = polyfit(*, h, n)Pm = polyval(P, *)R = sqrt(sum(h - Pm).2) %t = linspace(*(1), *(end), 12);Poly = polyval(P, t);figure(2)plot(*, 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 Appro*i

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

21、ewton-Cotes求积公式特别地，当取插值节点为时有 = 1 * GB3 两点公式梯形公式： = 2 * GB3 三点公式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 = 1 * GB2 试利用复化的梯形求积法求该组数据所在曲线与基准线轴在*围内所围成图形面积.画出数据散点图和图形的示意图. = 2 * GB2 试利用复化的Simpson求积法求该组数据所在曲线与基

22、准线轴在*围内所围成图形面积.画出数据散点图和图形的示意图.解答如下：程序：% *I = 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(*I) - 1; %n=10M = 60;* = linspace(*I(1), *I(end), M +1); % 1.1 T1 = 0; %for k = 2 : n + 1 T1 = T1 + (YI(k) + YI(k-1)*(*I(k) - *I(k

23、-1)/2;enddisp()T_Trape = T1 %figure(1) %set(gca,FontSize, 20)patch(*I(1), *I, *I(end), 0, YI, 0, 0.5 1 0.5) %hold onplot(*I, YI, o-, *I, 0*I, k, LineWidth, 3, markersize, 8)title()*label(*)ylabel(y)运行结果：复化的梯形求积法求得的图形面积：T_Trape = -11.6090图像：程序：函数M文件：function Ln = Lagrange_Fun_01(*, *I, YI)if isa(*, s

24、ym) = 1; n = length(*I) - 1; Ln = 0; Pn = sym(ones(n + 1, 1);for k = 1 : n + 1for i = 1 : n + 1if i = k Pn(k) = Pn(k)*(* - *I(i)/(*I(k) - *I(i);else Pn(k) = Pn(k);endend Ln = Ln + Pn(k)*YI(k);endelse n = length(*I) - 1; L = ones(n + 1, length(*); Ln = zeros(size(*);for k = 1 : n + 1for i = 1 : n + 1

25、if i = k L(k,:) = L(k,:).*(* - *I(i)./(*I(k) - *I(i);else L(k,:) = L(k,:);endend Ln = Ln + L(k, :).*YI(k);endend主程序：% 1.2 1.3%32S2 = 0;%Simpsonfor k = 1: 2: n - 1 S2 = S2 + (*I(k+2) - *I(k)*(YI(k) + 4*YI(k+1) + YI(k+2)/6; *(k, :) = linspace(*I(k), *I(k+2), 20); Y(k, :) = Lagrange_Fun_01(*(k,:), *I(k

26、 : k+2), YI(k : k+2);enddisp(Simpson)S_Simpson = S2 %Simpsonfigure(2)set(gca,FontSize, 20)hold onfor k = 1: 2 : n-1 patch(*I(k), *(k,:), *I(k +2), 0, Y(k, :), 0, 0.5 1 0.5) plot(*(k, :), Y(k, :), LineWidth, 3)endplot(*I, YI, o, *I, 0*I, k, LineWidth, 3, markersize, 8)title(Simpson)*label(*)ylabel(y)

27、运行结果：Simpson复化求积法求得的图形面积：S_Simpson = -11.9128图像：第八章微分方程组初值问题数值方法引言：微分方程组是科学研究和工程应用中最常用的数学模型之一。但是基绝大多数的常微分方程和常微分方程组得不到解析解和实际应用中往往只需要知道常微分方程组的解在*些点处的函数值这两点原因，绝大多数在实际中遇到的常微分方程和常微分方程组得不到解析解，我们可以采用下面将介绍的常微分方程组的初值问题的数值解法，就可以到达这一目的。Euler 方法最简单的数值解法1、Euler 法假设要求在点,，处初值问题4的解的近似值。首先对式4的两端积分，得对于该式的右边，如果被积函数用积分

28、下限处的函数值代替被积函数作积分，则有进而得到下式给出的递推算法Euler 方法2、Euler 方法的改良改良的Euler 方法Euler 方法的精度高，其原因在于：在推导Euler 方法时，我们是用待求解函数在一点处的变化率1代替在区间上的平均变化率；而在推导改良的Euler 方法时，我们是用待求解函数在两点处变化率的平均值2代替在区间上的平均变化率；显然，通常式2比式1更接近于在区间上的平均变化率。由此，人们受到启发：适当地选取区间上函数多个点处的变化率，用它们加权平均值代替在区间上的平均变化率，近似解的精度应更高。RungeKutta法： = 1 * GB3 二阶的RungeKutta法 = 2 * GB3 三阶的RungeKutta法 = 3 * GB3 四阶的RungeKutta法RungeKutta法是按选取区间上函数变化率数目的多少和截断误差的阶数来区分的一系列方法。例题：给定如下常微分方程Cauchy 问题：分别取参数，初值和步长：，。取，在区间上分别利用Euler 法，4阶R-K公式和matlab的ode45求问题1的数值解；程序：clcclear alla=1;b=1;F = (t, y) -a*y-t2*e*p(-0.5*t)*sin(t)+b;%matlabode45y0 = -1;Tf = 2*pi;Tspan = 0,