数值分析课设.doc

目录实验一 线性方程组数值解法1一、实验目的1二、实验基本原理1三、程序编写2四、实例分析3实验二 插值法和数据拟合3一、实验目的3二、实验基本原理3三、程序编写5四、实例分析5实验三 数值积分6一、实验目的6二、实验基本原理6三、程序编写7四、实例分析8实验四 常微分方程数值解8一、实验目的9二、实验基本原理9三、程序编写10四、实例分析10实验五 数值方法实际应用11一、实验目的11二、实验内容11三、实验基本原理11四、实例分析1217实验一 线性方程组数值解法一、实验目的1 了解LU分解法的基本原理和方法；2 通过实例掌握用MATLAB求线性方程组数值解的方法；3 编程实现LU分解二、实验基本原理2.1 LU分解法对于矩阵A，若存在一个单位下三角矩阵L和一个上三角U，使得（1.1） 即 （1.2）称上述分解为矩阵A的LU分解，或Doolittle分解，也称为直接三角分解。在式（1.2）中，利用矩阵L的第一行与矩阵U的各列相乘，可以得到矩阵U的第1行 （1.3）利用矩阵L的各行与矩阵U的第1列相乘，得到矩阵L的第1列 （1.4）假设已确定出矩阵U的第1行到第k-1行与矩阵L的第1列到第k-1列，现在来求矩阵U的第k行和L的第k列。利用式（1.2）中矩阵L的第k行与矩阵U的第列相乘，得到矩阵U的第k行 （1.5）利用矩阵L的第行与矩阵U的第k列相乘，得到矩阵L的第k列 （1.6） 显然，式（1.5）和式（1.6）对于都成立。若矩阵A有分解：A=LU，其中L为单位下三角阵，U为上三角阵，则称该分解为Doolittle分解，可以证明，当A的各阶顺序主子式均不为零时，Doolittle分解可以实现并且唯一。2.2.用矩阵的LU分解求解方程组由式（1.1），将方程组 改写为 则方程组求解可分成两部分完成。令 ，则方程组可改写成方程组和由上式得到 三、程序编写3.1矩阵的LU分解functionL,U,index=LU_Decom(A)n,m=size(A);if n=m error(The rows and columns of matrix A must be equal!); return;endL=eye(n);U=zeros(n);index=1;for k=1:n for j=k:n z=0; for q=1:k-1 z=z+L(k,q)*U(q,j); end U(k,j)=A(k,j)-z; end if abs(U(k,k) A=2 5 -6;4 13 -19;-6 -3 -6; L,U,index=LU_Decom(A)得到L = 1 0 0 2 1 0 -3 4 1U = 2 5 -6 0 3 -7 0 0 4index = 1index = 1表明计算成功，得到相应的分解矩阵。用矩阵的LU分解求解方程组在Matlab中输入： A=2 5 -6;4 13 -19;-6 -3 -6;b=10 19 -30; x,y=bxzylu(A,b)得到y = 10 -1 4x = 3 2 1实验二 插值法和数据拟合一、 实验目的1 了解lagrange插值法的基本原理和方法；2 通过实例掌握用MATLAB求插值的方法；3 编程实现lagrange插值二、 实验基本原理已知n+1个互异的差值节点的函数值，构造次数不超过n的多项式。令 （2.1）式中，是次插值基函数，并且满足 （2.2）由式（2.1）和式（2.2）可知，所得到的满足插值条件 。对于固定的，由式（2.2）可知，以为零点，并满足，由于是次多项式，经推导得到 （2.3）将式（2.3）代入式（2.1）得到 （2.4）称式（2.4）为次Lgrange插值公式。引进基函数则有 基函数表示 还有一种表示式 n+1(x) 表示式 其中。从而就得到了在区间a,b上的关于函数 f(x) 的n次多项式近似计算式： 三、 程序编写1.Lagrange插值（lagrange.m）function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for k=1:m z=x(k); s=0.0; for j=1:n p=1.0; for i=1:n if i=j p=p*(z-x0(i)/(x0(j)-x0(i); end end s=p*y0(j)+s; end y(k)=s;end2.主程序Lagrange2调用子程序x0=1 23 4 6 8 10 12 14 16; y0=4.00 6.41 8.01 8.79 9.53 9.86 10.33 10.42 10.53 10.61;x=1:0.1:16; y=lagrange(x0,y0,x); x1=5;x2=16.4;y1=lagrange(x0,y0,x1)y2=lagrange(x0,y0,x2)plot(x0,y0,.k,markersize,20)hold onplot(x,y,-r,markersize,30)hold on plot(x1,y1,*b,markersize,8)hold on plot(x2,y2,*b,markersize,8)legend(原数值点,Lagrange曲线,2)四、 实例分析某化学反应中,在有限个时刻t(min)，测得生成物质量浓度 y(10-3g/cm3)的如下数据ti12346810121416yi4.006.418.018.799.539.8610.3310.4210.5310.61那么在时刻t=5min，t=16.4min时的浓度是多少？解：编写程序进行计算(1). 先编写Lagrange插值算法子程序(2). 再编写主程序Lagrange2调用子程序。(3). 计算在t=5，t=16.4的浓度值，并画出曲线图。计算结果： y1=lagrange(x0,y0,x1)y1 = 9.2300 y2=lagrange(x0,y0,x2)y2 = 6.5872结果表明：t=5时，浓度值在Lagrange曲线上；t=16.4时，浓度值不在曲线上。说明内插比外插的精确度高。实验三 数值积分一、 实验目的1 了解数值积分的基本原理和方法；2 掌握用MATLAB求积分的方法；3 通过实例学习如何用几种方法求积分。4 了解书上介绍的几种数值积分的不同原理和方法；二、实验基本原理2.1 复化Simpson求积公式的原理将区间a,b分为n等分，在每个子区间上采用辛普森公式： ，若记，则得：。则可以记得：，称为复化辛普森求积公式。其余项为：，此外，由于中求积系数均为正数，故知复化辛普森公式计算稳定。2.2 复化梯形求积公式的原理将区间a,b分为n等分，分点在每个子区间上采用梯形公式计算，则得：，记：，称为复化梯形公式。其余项为：。三、 程序编写3.1 复化Simpson求积公式function I=S_quad(x,y)n=length(x);m=length(y);if n=m erroe(The lengths of X and Y must be equal! ); return;endif rem(n-1,2)=0 %如果n-1不能被2整除，则调用复化梯形公式 I=T_quad(x,y); return;endN=(n-1)/2;h=(x(n)-x(1)/N;a=zeros(1,n);for k=1:N a(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);3.2 复化梯形求积公式function I=T_quad(x,y)n=length(x);m=length(y);if n=m erroe(The lengths of X and Y must be equal! ); return;endh=(x(n)-x(1)/(n-1);a=1 2*ones(1,n-2) 1;I=h/2*sum(a.*y);四、 实例分析分别利用复化梯形公式、复化simpson 公式计算积分，在积分区间中点与点之间的间隔取为0.1。解：（1）用复化梯形公式求解在Matlab中输入： x=-1:0.1:1; y=exp(-x.2); I=T_quad(x,y)得到结果：I = 1.4924（2）用复化simpson 公式求解在Matlab中输入： x=-1:0.1:1; y=exp(-x.2); I=S_quad(x,y)得到结果：I = 1.4936实验四 常微分方程数值解一、 实验目的1了解matlab中提供的求解常微分方程解的常见方法；2 了解Runge-Kutta法的收敛性与稳定性；3 学会用Matlab编程实现Runge-Kutta法解常微分方程；4 熟记Runge-Kutta法的各阶表达形式；二、实验基本原理2.1 改进Euler公式在计算中，先用Euler公式求得一个初步的近似值（称为预测值），再用梯形公式将它校正一次，得，这称为预测校正公式，即改进Euler方法，其迭代格式为预测：校正：2.2 龙格-库塔(Runge-Kutta)方法龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：,使用差分概念。推出（近似等于，极限为）另外根据Lagrange微分中值定理，存在,使得其中，称为在上的平均斜率。利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出常用的四阶龙格库塔公式：4阶经典Runge-Kutta格式 三、 程序编写3.1 改进Euler公式functionx,y=Euler_r(ydot_fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h; ybar=y(n)+h*feval(ydot_fun,x(n),y(n); y(n+1)=y(n)+h/2*(feval(ydot_fun,x(n),y(n). +feval(ydot_fun,x(n+1),ybar);end3.2 四阶龙格-库塔(Runge-Kutta)方法functionx,y=Runge_Kutta4(ydot_fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:N x(n+1)=x(n)+h; k1=h*feval(ydot_fun,x(n),y(n); k2=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k1); k3=h*feval(ydot_fun,x(n)+1/2*h,y(n)+1/2*k2); k4=h*feval(ydot_fun,x(n)+h,y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end四、 实例分析分别用改进Euler方法和标准四阶Runge-Kutta方法求初值问题的解在处的近似值，步长分别取 和。解：（1）改进Euler方法取 h=0.02，N=5在Matlab中输入：ydot_fun=inline(x+y,x,y); x,y=Euler_r(ydot_fun,0,1,0.02,5)得到结果：x = 0 0.0200 0.0400 0.0600 0.0800 0.1000y = 1.0000 1.0204 1.0416 1.0637 1.0866 1.1103（2）标准四阶Runge-Kutta方法在Matlab中输入： ydot_fun=inline(x+y,x,y); x,y=Runge_Kutta4(ydot_fun,0,1,0.1,1)得到结果：x = 0 0.1000y = 1.0000 1.1103四阶Runge-Kutta方法中取 h=0.1达到改进Euler方法中 h=0.02的计算结果（实际上更精确，只是显示结果相同），这表明，四阶Runge-Kutta方法有很好的计算效果。实验五 数值方法实际应用一、实验目的1 了解最小二乘拟合的基本原理和方法；2 掌握用MATLAB作最小二乘多项式拟合和曲线拟合的方法；3 通过实例学习如何用拟合方法解决实际问题，注意与插值方法的区别。4 了解各种参数辨识的原理和方法；5 通过范例展现由机理分析确定模型结构，拟合方法辨识参数，误差分析等求解实际问题的过程；二、实验内容1用MATLAB中的函数作一元函数的多项式拟合与曲线拟合；2用MATLAB中的函数作二元函数的最小二乘拟合，作出误差图；3针对预测和确定参数的实际问题，建立数学模型，并求解。三、实验基本原理最小二乘法的思想是：设有N个数据点（xk,yk），并给定M个线性独立函数fj(x)。为求M个稀疏，使用由线性组合形成的函数f(x)，表示为：求解最小误差平方和，表示为：为求解E的最小值，每个偏导数必须为零（即），这样可得到如下方程组：交换上式中的求和顺序，可得一个M*M线性方程组，未知数是系数cj，称为正规方程组： 将上面的正规方程组表示成矩阵形式可减少不必要的计算量。方程组的系数矩阵是一个对称矩阵，并且是正定的，可以唯一地确定系数cj，根据这些特点，构造如下矩阵： 乘积FF是一个M*M矩阵，方程组转化为求解线性方程组：（求解系数矩阵C）四、实例分析4.1 问题重述根据实验内容和步骤，完成以下具体实验 旧车价格预测某年美国旧车价格的调查资料如下表，其中xi表示轿车的使用年数，yi表示相应的平均价格。试分析用什么形式的曲线来拟合上述的数据，并预测使用4.5年后轿车的平均价格大致为多少？xi12345678910yi2615194314941087765538484290226204多项式函数拟合：a,s=polyfit(xdata,ydata,n)其中n表示多项式的最高阶数，xdata，ydata为将要拟合的数据，它是用数组的方式输入输出参数a为拟合多项式的系数,s为误差。另：多项式在x处的值y可用下面程序计算y=polyval(a,x)，对上面给出的数据做多项式拟合，可取不同的n观测此时的误差，看取什么样的n较好。4.2求解过程首先画出数据的散点图。输入 x=1 2 3 4 5 6 7 8 9 10;y=2615 1943 1494 1087 765 538 484 290 226 204; plot(x,y,+)进行多项式拟合：编写Matlab程序function p=mafit(x,y,m)format short;A=zeros(m+1,m+1);for i=0:m for j=0:m A(i+1,j+1)=sum(x.(i+j); end b(i+1)=sum(x.i.*y);enda=Ab;p=fliplr(a);运行结果为： x=1 2 3 4 5 6 7 8 9 10;y=2615 1943 1494 1087 765 538 484 290 226 204; p=mafit(x,y,2)p = 1.0e+003 * 0.0361 -0.6508 3.1523 p=mafit(x,y,3)p = 1.0e+003 * -0.0027 0.0800 -0.8528 3.3801 p=mafit(x,y,4)p = 1.0e+003 * 0.0001 -0.0054 0.1002 -0.9082 3.4232 p=mafit(x,y,5)p = 1.0e+003 * 0.0000 -0.0012 0.0078 0.0410 -0.7954 3.3549 p=mafit(x,y,6)p = 1.0e+003 * 0.0001 -0.0033 0.0411 -0.2577 0.8811 -2.0266 3.9789 p=mafit(x,y,7)p = 1.0e+003 * Columns 1 through 7 -0.0000 0.0016 -0.0261 0.2261 -1.0888 2.9137 -4.4740 Column 8 5.0632 p=mafit(x,y,8)p = 1.0e+003 * Columns 1 through 7 -0.0000 0.0013 -0.0228 0.2171 -1.2001 3.8959 -7.0902 Columns 8 through 9 5.8814 0.9322整理得到拟合多项式：2阶拟合多项式：3阶拟合多项式：4阶拟合多项式：5阶拟合多项式：6阶拟合多项式：7阶拟合多项式：8阶拟合多项式：为检验每个拟合多项式的效果，编写程序lspoly2.m：function Y=lspoly2(m)x=1:1:10;y=2615 1943 1494 1087 765 538 484 290 226 204;C=(mafit(x,y,m); n=length(x);X=zeros(n,m+1);for k=1:m+1 X(:,k)=x.(k-1);endY=(X*C);运行结果为： Y=lspoly2