数值分析课程设计.doc

摘要实验一 拉格朗日插及数值求解1.1实验目的了解Lagranger差值的基本原理和方法通过实例掌握用MATLAB求插值的方法根据实际计算理论，利用Lagranger插值多项式计算1.2实验原理设已知，.,及=f()(i=0,1,.,n),为不超过n次多项式且满足（i=0,1,.n）.易知其中，均为n次多项式,再由（ji）为n次多项式的n个根知.最后，由,i=0,1,.,n.总之，=，=式为n阶Lagrange插值公式，其中，（i=0,1,.n）称为n阶Lagrange插值的基函数。1.3实验内容function y = lagranger(x0,y0,x);%UNTITLED Summary of this function goes here% Detailed explanation goes heren=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n li=1.0; for j=1:n if j=k li=li*(z-x0(j)/(x0(k)-x0(j); end end s=li*y0(k)+s; end y(i)=s; end1.4实验案例及结果分析（1）输入：x0=4,5,6;y0=10,5.25,1;x=5;y=lagranger(x0,y0,x)（2）输入：X0=1,4,8;y0=6,3.2,4;x=4;y=lagranger(x0,y0,x)实验二 LU分解法解线性方程组2.1实验目的1.了解LU分解法解线性方程组的基本原理；2.熟悉计算方法的技巧和过程，能用LU分解法解实际问题；3.用matlab实现LU分解。2.2实验原理1.若一个线性方程组系数矩阵为n阶方阵A且各阶顺序主子式均不为0则A的LU分解存在且唯一。将高斯消去法改写为紧凑形式，可以直接从矩阵A的元素得到计算L，U元素的递推公式，而不需任何中间步骤，这就是所谓直接三角分解法。一旦实现了矩阵A的LU分解，那么求解Ax=b的问题就等价于求解两个三角形方程组：Ly=b，求y；Ux=y，求x。2.在满足1的条件下课推导得出以下公式（1） （2） 3.公式（1）用于求解矩阵L、U，公式（2）用于会带求解y、x。从公式中可以看出：L对角线上元素为1，U第一行与A第一行相同。4LU分解的具体过程和顺序如下：（1）第一步分解：（2）第二步分解： （3）第三步分解： （n）第n步分解：依次计算：、.，. 2.3实验内容编写一个M文件function L,U,x=Lu_x(A,b)n,m=size(A);if n=merror(The rows and columns of matrix A must be equal!);return;endfor ii=1:nfor i=1:iifor j=1:iiAA(i,j)=A(i,j);endendif (det(AA)=0)error(The matrix can not be divided by LU!)return;endendA n,n=size(A); L=zeros(n,n); U=zeros(n,n); for i=1:n L(i,i)=1; end for k=1:n for j=k:n U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j); end for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)/U(k,k); end end y(1)=d(1);for i=2:n for j=1:i-1d(i)=d(i)-L(i,j)*y(j);endy(i)=d(i);endx(n)=y(n)/U(n,n);for i=(n-1):-1:1for j=n:-1:i+1y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);end2.4实验案例及结果分析在MATLAB命令窗体输入：A=10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2;b=8,5.900001,5,1; L,U,x=Lu_x(A,b)得到结果如下：实验三 龙贝格求积公式求数值积分3.1实验目的1. 熟练掌握龙贝格求积的基本思路和步骤；2. 培养编程与上机调试能力；3.利用龙贝格求积方法求解积分。3.2实验原理（1）置n=1，精度要求，（2）计算（3）置，并计算（4）置m=n，n=2n，k=1。（5）计算。（6）若m=1，转（7）；否则，置，k=k+1，转（5）。（7）若，则停止计算（输出），否则转（3）。3.3实验内容function R = romberg(f, a, b, e)% 参数介绍:% f - 被积函数f(x)% a - x的左区间.% b - x的右区间.% e - 误差限.% 结果:% R - 返回Romberg 表. n = 1; %区间二分次数while 1 %在此仅代表多次二分，在后面判断循环终止 R = zeros(n + 1, n + 1);%生成(n+1)*(n+1)的0矩阵 R(0+1, 0+1) = (b - a) / 2 * (feval(f, a) + feval(f, b); % 初始值(2点梯形公式). for i = 1 : n % 按照公式计算Romberg 表的第一列. h = (b - a) / 2i; s = 0; for k = 1 : 2(i-1) s = s + feval(f, a + (2*k - 1)*h); end R(i+1, 0+1) = R(i-1+1, 0+1)/2 + h*s; end for j = 1 : n % 计算Romberg 表的其他列. fac = 1 / (4j - 1); for m = j : n R(m+1, j+1) = R(m+1, j-1+1) + fac*(R(m+1, j-1+1) - R(m-1+1, j-1+1); end end if abs(R(n,n) - R(n+1,n+1) e %当精度满足设定要求时退出 break; end n = n + 1; %未达到指定精度继续二分endt=R(i,j)3.4实验案例及结果分析用Romberg求积法计算积分，取精度要求=10-5。在Matlab命令窗口输入：fun=inline(5./(3+x),x) ;romberg(fun,-1,1,1e-6)输出结果如下：由输出结果可知最终积分为3.4657。实验四 Runge-Kutta方法求常微分方程数值解4.1实验目的1.熟悉Runge-Kutta常微分方程初值问题的基本原理2.了解Runge-Kutta常微分方程初值问题的计算流程3.能编程实现Runge-Kutta常微分方程初值问题4.2实验原理在欧拉法的基础上增加了计算一个右函数f的值，可望p=2。若要使得到的公式阶数p更大，就必须包含更多的f值。实际上从与方程等价的积分形式即若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必须要增加求积节点，为此可将上述公式的右端项表示为一般来说，点数r越多，精度越高，上式右端相当于增量函数，为得到便于计算的显示方法，可类似改进的欧拉法，将公式表示为这里均为常数，上式称为显示龙格库塔（Ruuge-Kutta）法，简称R-K方法。当r=1时，就是欧拉法，当r=2时，改进的欧拉法就是其中的一种。依次类推，如果在区间内多预估几个点上的斜率值并用他们的加权平均数作为平均斜率的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格库塔公式，也就是在工程中应用广泛的经典龙格库塔算法： 4.3实验内容四阶龙格库塔法的计算公式为：四阶龙格库塔公式的Matlab程序代码：function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2); K3 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K2*h/2); K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*K3); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;end函数运行时需要调用下列函数： function fv=Funval(f, varvec, varval) var= findsym(f); if length(var) %多项式插值 t=205,430,677,945,1233,1542,1872,2224; c=100,200,300,400,500,600,700,800; c1=0:10:1000; t1=interp1(c,t,c1,nearest); t2=interp1(c,t,c1,linear); t3=interp1(c,t,c1,spline); t4=interp1(c,t,c1,cubic); plot(c,t,*,c1,t1,:b,c1,t2,-r,c1,t3,-g,c1,t4,.-r) legend(原始数据,最近点插值,线性插值,样条插值,立方插值)输出结果如下： %多项式拟合 t=205,430,677,945,1233,1542,1872,2224; c=100,200,300,400,500,600,700,800; n=8; %拟合度 p=polyfit(c,t,n)Warning: Polynomial is not unique; degree = number of data points. In polyfit at 72p = 1.0e-006 * 0.0000 -0.0000 0.0000 -0.0048 0.8068 0 0 0 0 c1=linspace(0,10,1000);%x轴范围 t1=polyval(p,c1);%计算在c1数据点的多项式值 plot(c,t,*,c,t,c1,t1,:b) xlabel(c); ylabel(t=t(c); title(八次曲线拟合)图2 八次多项式拟合图通过实验结果可知:(1)播放时间随计数器读数增大而增大；(2)与计数器读数具有近似线性关系；总结通过本次实验，我可以利用matlab的一些基本功能对问题进行求解，初步利用Matlab进行实验操作，通过此次实验也使我了解到MATLAB数学软件功能的强大，学会了用MATLAB进行编程，解决数值分析中的实际问题，同时在编程过程中，更加深入的了解了LU分解法，拉格朗日插值，龙贝格求积，龙格库塔的原理，对算法基本原理认识深刻，进一步加深相关算法的效用和使用范围的全面了解。没有像书中那样进行枯燥的定理论述，而是重点介绍了方法基本思想以及在MATLAB平台上的使用。课堂气氛更加活跃，使我们在不知不觉中掌握了大部分知识。还学会了设计算法不单单是设计出来，一个好的算法就要保证算法具有收敛性和较高的收敛速度，保证算法具有数值的稳定性，小