数值分析的matlab实现
MATLAB数值实验一(数据的插值运算及其应用完整版)
佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 插值法与数据拟合 专业班级 机械工程 姓 名 余红杰 学 号 10 指导教师 陈剑 成 绩 日 期 月 日一、实验目的1、学会Lagrange 插值、牛顿插值和三次样条插值等基本插值方法；2、讨论插值的Runge 现象3、学会Matlab 提供的插值函数的使用方法，会用这些函数解决实际问题。
二、实验原理1、拉格朗日插值多项式2、牛顿插值多项式3、三次样条插值 三、实验步骤1、用MATLAB 编写独立的拉格朗日插值多项式函数2、用MATLAB 编写独立的牛顿插值多项式函数3、用MATLAB 编写独立的三次样条函数（边界条件为第一、二种情形）4、已知函数在下列各点的值为：根据步骤1,2,3编好的程序，试分别用4次拉格朗日多项式4()L x 、牛顿插值多项式4()P x 以及三次样条函数()S x (自然边界条件)对数据进行插值，并用图给出 ｛(,),0.20.08,0,1,2,,10i i i x y x i i =+=｝，4()L x 、4()P x 和()S x 。
5、在区间[-1,1]上分别取10,20n =用两组等距节点对龙格函数21(),(11)125f x x x=-≤≤+作多项式插值，对不同n 值，分别画出插值函数及()f x 的图形。
6、下列数据点的插值可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9个点作8次多项式插值8()L x 。
（2）用三次样条（第一边界条件）程序求()S x 。
7、对于给函数21()125f x x =+在区间[-1,1]上取10.2(0,1,,10)i x i i =-+=，试求3次曲线拟合，试画出拟合曲线并打印出方程，与第5题的结果比较。
四、实验过程与结果：1、Lagrange 插值多项式源代码：function ya=lag(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 ya=0; mu=1; %初始化%循环方式求L 系数，并求和： for i = 1:length(y) for j = 1:length(x) if i ~= jmu = mu * (xa - x(j) ) / ( x(i) - x(j) ); else continue end endya = ya + y(i) * mu ; mu = 1; end2、Newton 源代码：function ya = newton(x,y,xa) %x 所有已知插值点 %y 插值点对应函数值 %xa 所求点，自变量 %ya 所求点插值估计量 %建立系数零矩阵D 及初始化：D = zeros(length(x)-1);ya = y(1);xi = 1;%求出矩阵D，该矩阵第一行为牛顿插值多项式系数：for i=1:(length(x)-1)D(i,1) = (y(i+1) -y(i))/(x(i+1) -x(i));endfor j=2:(length(x)-1)for i=1:(length(x)-j)D(i,j) = (D(i+1,j-1) - D(i,j-1)) / (x(i+j) - x(i)); endend%xi为单个多项式（x-x(1))(x-x(2))...的值for i=1:(length(x)-1)for j=1:ixi = xi*(xa - x(j));endya = ya + D(1,i)*xi;xi = 1;end3、三次样条插值多项式（1）（第一边界条件）源代码：function y=yt1(x0,y0,f_0,f_n,x) _____________(1)%第一类边界条件下三次样条插值；%xi 所求点；%yi 所求点函数值；%x 已知插值点；%y 已知插值点函数值；%f_0左端点一次导数值；%f_n右端点一次导数值；n = length(x0);z = length(y0);h = zeros(n-1,1);k=zeros(n-2,1);l=zeros(n-2,1);S=2*eye(n);for i=1:n-1h(i)= x0(i+1)-x0(i);endfor i=1:n-2k(i)= h(i+1)/(h(i+1)+h(i));l(i)= 1-k(i);end%对于第一种边界条件：k = [1;k]; _______________________(2)l = [l;1]; _______________________(3)%构建系数矩阵S：for i = 1:n-1S(i,i+1) = k(i);S(i+1,i) = l(i);end%建立均差表：F=zeros(n-1,2);for i = 1:n-1F(i,1) = (y0(i+1)-y0(i))/(x0(i+1)-x0(i));endD = zeros(n-2,1);for i = 1:n-2F(i,2) = (F(i+1,1)-F(i,1))/(x0(i+2)-x0(i));D(i,1) = 6 * F(i,2);end%构建函数D：d0 = 6*(F(1,2)-f_0)/h(1); ___________(4)dn = 6*(f_n-F(n-1,2))/h(n-1); ___________(5)D = [d0;D;dn]; ______________(6)m= S\D;%寻找x所在位置，并求出对应插值：for i = 1:length(x)for j = 1:n-1if (x(i)<=x0(j+1))&(x(i)>=x0(j))y(i) =( m(j)*(x0(j+1)-x(i))^3)/(6*h(j))+...(m(j+1)*(x(i)-x0(j))^3)/(6*h(j))+...(y0(j)-(m(j)*h(j)^2)/6)*(x0(j+1)-x(i))/h(j)+... (y0(j+1)-(m(j+1)*h(j)^2)/6)*(x(i)-x0(j))/h(j) ; break;else continue;endendend（2）（自然边界条件）源代码：仅仅需要对上面部分标注的位置做如下修改:__(1):function y=yt2(x0,y0,x)__(2):k=[0;k]__(3):l=[l;0]__(4)+(5):删除—(6):D=[0:D:0]4、——————————————PS:另建了一个f方程文件，后面有一题也有用到。
有限差分法matlab程序
有限差分法matlab程序
有限差分法是一种常见的数值分析方法，它可以用来求解微分方程的解，称为数值解。
有限差分法的实现可以使用MATLAB，MATLAB 的强大的编程能力可以帮助我们更方便的实现有限差分法。
首先，我们需要了解什么是有限差分法。
有限差分法是一种数值计算的方法，它用一组等间距的数组来拟合一个函数，然后利用这些数组求出函数的值。
有限差分法可以用来解决复杂的微分方程，这是由于它能够将一个复杂的微分方程分解成一系列简单的微分方程，每个方程都可以用有限差分法来解决。
其次，使用MATLAB来实现有限差分法需要了解MATLAB中有限差分法的必要组件。
MATLAB中有限差分法的组件主要有：矩阵操作、函数计算、微分方程求解等。
首先，我们使用矩阵操作来实现函数的拟合，然后使用函数计算得到函数的值。
最后，我们使用微分方程求解器来解决复杂的微分方程。
最后，实现有限差分法在MATLAB中需要编写MATLAB程序。
MATLAB 程序首先需要编写矩阵操作程序，使用矩阵操作程序拟合函数；然后编写函数计算程序，计算函数的值；最后，编写微分方程求解程序，利用微分方程求解器来解决复杂的微分方程。
有限差分法是一种重要的数值分析方法，在很多科学中都有着广泛的应用，MATLAB可以让我们更方便的实现有限差分法，在微分方程的求解中发挥其独特的作用。
本文简要介绍了如何使用MATLAB来实现有限差分法，希望能为大家提供一些帮助。
第十章MATLAB的数值分析
• 第一个问题可归结为“已知函数在x0,x1,
– …,xn处的值，求函数在区间[x0,xn]内其它点处的值”，这 种问题适宜用插值方法解决。 – 插值问题可描述为：已知函数在x0,x1,…,xn处的值 y0,y1,…,yn，求函数p(x)，使p(xi) = yi。
• 但对第二个问题不宜用插值方法,因为600米已超出所 给数据范围，用插值函数外推插值区间外的数据会 产生较大的误差。
– Q1=prctile(w,25); – Q3=prctile(w,75); – prctile( )函数实现计算样本的百分位数功能
分布形态的测定
• 只用集中趋势和离中趋势来表示所有数据，难免不 够准确。分析总体次数的分布形态有助于识别整个 总体的数量特征。总体的分布形态可以从两个角度 考虑，一是分布的对称程度，另一个是分布的高低。 前者的测定参数称为偏度或偏斜度，后者的测定参 数称为峰度。 • 峰度是掌握分布形态的另一指标，它能描述分布的 平缓或陡峭程度。如果峰度数值等于零，说明分布 为正态；若峰度数值大于零，说明分布呈陡峭状态； 若峰度数值小于零，说明分布形态趋于平缓。
– 解决第二个问题的常用方法是，根据地面到井下 500 处的 数据求出瓦斯浓度与地面到井下距离x之间的近似函数关 系f(x), 由f(x)求井下600米处的瓦斯浓度。
• 插值函数过已知点，拟合函数不一定过已知点。通 常, 插值主要用于求函数值，而拟合的主要目的是求 函数关系。当然，某些问题既可以用插值也可以用 拟合。
插值方法-概述
• 为什么需要插值？
(1) 函数关系y=f(x)没有明确的表达式
(2) y=f(x)表达式复杂，不便于研究和使用
-20 -15
沉陷量/mm 下沉方向为"+"
数值分析(hilbert矩阵)病态线性方程组的求解matlab程序
（Hilbert 矩阵）病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n Hh ，11ij h i j ，i ，j = 1,2，…，n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ，取(1,1,,1)nx K ?，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，比较结果。
3.结合计算结果，试讨论病态线性方程组的求解。
第1小题：condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出，1）在n小于等于13之前，图像近似直线log（cond（H））~1.519n-1.8332）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势第2小题：solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);%消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多，迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多，迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多，迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多，迭代可能不收敛');break;endend。
牛顿插值法matlab程序例题
牛顿插值法是一种常用的数值分析方法，用于构造一个多项式函数，以便在给定的数据点上进行插值。
这个主题在数学和工程领域中有着广泛的应用，特别是在数据拟合和函数逼近方面。
牛顿插值法的核心思想是通过不断地添加新的数据点来构造一个多项式，并利用已知数据点来确定多项式的系数，从而实现对未知数据点的插值预测。
在Matlab中，实现牛顿插值法并不困难，我们可以利用已有的函数和工具来简化计算过程。
下面，我们将通过一个具体的例题来讲解如何使用Matlab编写牛顿插值法的程序，并分析其结果。
我们需要明确牛顿插值法的数学原理。
给定n个互不相同的节点\(x_0, x_1, ... , x_n\)，以及在这些节点上的函数值\(f(x_0), f(x_1), ... , f(x_n)\)，我们希望构造一个n次插值多项式p(x)，满足p(x_i) = f(x_i)，i=0,1,...,n。
牛顿插值多项式的一般形式为：\[p(x) = a_0 + a_1(x - x_0) + a_2(x - x_0)(x - x_1) + ... + a_n(x -x_0)(x - x_1)...(x - x_{n-1})\]其中，\[a_i\]表示插值多项式的系数。
通过牛顿插值法的迭代过程，可以逐步求解出这些系数，进而得到插值多项式的表达式。
接下来，我们将以一个具体的例题来演示如何在Matlab中实现牛顿插值法。
假设我们有如下的数据点和函数值：\(x = [1, 2, 3, 4]\)\(f(x) = [1, 4, 9, 16]\)我们希望利用这些数据点来构造一个插值多项式，并在给定的区间上进行插值计算。
在Matlab中，可以通过interp1函数来进行插值计算，该函数支持多种插值方法，包括牛顿插值法。
下面是一个简单的Matlab程序示例：```matlabx = [1, 2, 3, 4];y = [1, 4, 9, 16];xi = 2.5;yi = interp1(x, y, xi, 'spline');disp(['在x=',num2str(xi),'处的插值结果为：',num2str(yi)]);```在这段代码中，我们首先定义了给定的数据点x和对应的函数值y，然后利用interp1函数对x=2.5处的插值结果进行计算。
数值分析的matlab实现
第2章牛顿插值法实现参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55.求牛顿插值多项式和差商的MA TLAB 主程序：function[A,C,L,wcgs,Cw]=newpoly(X,Y)n=length(X);A=zeros(n,n);A(:,1) =Y';s=0.0;p=1.0;q=1.0;c1=1.0;for j=2:nfor i=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));endb=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1;endC=A(n,n);b=poly(X(n));q1=conv(q1,b);for k=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);endL(k,:)=poly2sym(C);Q=poly2sym(q1);syms Mwcgs=M*Q/c1;Cw=q1/c1;（1）保存名为newpoly.m 的M 文件（2）输入MA TLAB 程序>> X=[242,243,249,250];>> Y=[13.681,13.526,13.098,13.095];>> [A,C,L,wcgs,Cw]=newpoly(X,Y)输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其)()()1(ξ+n n f x R 的系数向量Cw 。
A =13.6810 0 0 013.5260 -0.1550 0 013.0980 -0.0713 0.0120 013.0950 -0.0030 0.0098 -0.0003C =1.0e+003 *-0.0000 0.0002 -0.0551 4.7634L =- (23*x^3)/84000 + (2981*x^2)/14000 - (7757472138947345*x)/140737488355328 + 5237382665812919/1099511627776wcgs =(M*(x^4 - 984*x^3 + 363071*x^2 - 59535444*x + 3660673500))/24Cw =1.0e+008 *0.0000 -0.0000 0.0002 -0.0248 1.5253输入MATLAB程序>> x=244;>> y=- (23*x^3)/84000 + (2981*x^2)/14000 - (7757472138947345*x)/140737488355328 + 5237382665812919/1099511627776y =13.3976输入MATLAB程序>> x=[244,245,246,247,248];>> y=- (23.*x.^3)./84000 + (2981.*x.^2)./14000 - (7757472138947345.*x)./140737488355328 + 5237382665812919./1099511627776y =13.3976 13.2943 13.2143 13.1560 13.1178第4章 高斯-勒让德积分公式实现用高斯-勒让德积分公式计算dx e I x ⎰--=112221π，取代数精度为3和5，再根据截断误差公式写出误差公式，并将计算结果与精确值进行比较。
数值分析试验幂法与反幂法matlab教程文件
数值分析试验幂法与反幂法m a t l a b一、问题的描述及算法设计（一）问题的描述我所要做的课题是：对称矩阵的条件数的求解设计1、求矩阵A 的二条件数问题 A=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----210121012 2、设计内容：1）采用幂法求出A 的. 2）采用反幂法求出A 的.3）计算A 的条件数 ⅡA Ⅱ2* ⅡA -1Ⅱ2=cond2（A ）=/.（精度要求为10-6）3、设计要求1）求出ⅡA Ⅱ2。
2）并进行一定的理论分析。
（二）算法设计1、幂法算法（1）取初始向量u )0(（例如取u )0(=(1,1,…1)T ）,置精度要求ε，置k=1.（2）计算v )(k =Au )1(-k ,m k =max(v )(k ), u )(k = v )(k / m k（3）若| m k = m 1-k |<ε，则停止计算（m k 作为绝对值最大特征值1λ，u )(k 作为相应的特征向量）否则置k=k+1，转（2）2、反幂法算法（1）取初始向量u )0(（例如取u )0(=(1,1,…1)T ）,置精度要求ε，置k=1.（2）对A 作LU 分解，即A=LU（3）解线性方程组 Ly )(k =u )1(-k ,Uv )(k =y )(k（4）计算m k =max(v )(k ), u )(k = v )(k / m k（5）若|mk =m1-k|<ε，则停止计算（1/mk作为绝对值最小特征值nλ，u)(k作为相应的特征向量）;否则置k=k+1，转（3）.二、算法的流程图（一）幂法算法的流程图（二）反幂法算法的流程图三、算法的理论依据及其推导（一）幂法算法的理论依据及推导幂法是用来确定矩阵的主特征值的一种迭代方法，也即，绝对值最大的特征值。
稍微修改该方法，也可以用来确定其他特征值。
幂法的一个很有用的特性是它不仅可以生成特征值，而且可以生成相应的特征向量。
实际上，幂法经常用来求通过其他方法确定的特征值的特征向量。
matlab实现数值分析插值及积分
数值分析学院：计算机专业：计算机科学与技术班级：xxx学号：xxx姓名：xxx指导教师：xxx数值分析摘要：数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。
在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。
学习数值分析这门课程可以让我们学到很多的数学建模方法。
分别运用matlab数学软件编程来解决插值问题和数值积分问题。
题目中的要求是计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。
编程求解出来的结果为：P4(x)=x4+1。
其中Aitken插值计算的结果图如下：对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。
编程求解出来的结果为：0.6932其中复化梯形公式计算的结果图如下：问题重述问题一：已知列表函数表格 1x 0 1 2 3 4 f (x )121782257分别用拉格朗日，牛顿，埃特金插值方法计算P 4(x )。
问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分∫1x 21dx ，使精度小于5×10−5。
问题解决问题一：插值方法对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。
一、拉格朗日插值法：拉格朗日插值多项式如下：首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数)(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子)())(()(110n i i x x x x x x x x ----+- 。
又因)(x l i 是一个次数不超过n 的多项式，所以只可能相差一个常数因子，固)(x l i 可表示成：)())(()()(110n i i i x x x x x x x x A x l ----=+-利用1)(=i i x l 得：)())(()(1110n i i i i i i x x x x x x x x A ----=+-于是),,2,1,0()())(()()())(()()(110110n i x x x x x x x x x x x x x x x x x l n i i i i i i n i i i =--------=+-+-因此满足i i n y x L =)( ),2,1,0(n i =的插值多项式可表示为：∑==nj j j n x l y x L 0)()(从而n 次拉格朗日插值多项式为：),,2,1,0()()(0n i x l y x L nj i j j i n ==∑=matlab 编程：编程思想：主要从上述朗格朗日公式入手：依靠循环，运用poly （）函数和conv （）函数表示拉格朗日公式，其中的poly （i ）函数表示以i 作为根的多项式的系数，例如poly （1）表示x-1的系数，输出为1 -1，而poly （poly （1））表示（x-1）*（x-1）=x^2-2*x+1的系数，输出为1 -2 1；而conv （）表示多项式系数乘积的结果，例如conv （poly （1），poly （1））输出为1 -2 1；所以程序最后结果为x^n+x^n-1+……+x^2+x+1（n 的值据结果的长度为准）的对应系数。
数值分析matlab上机实验报告
数值分析matlab上机实验报告matlab软件实验报告数学上机课实验报告matlab实验报告总结数值分析试卷篇一：《MATLAB与数值分析》第一次上机实验报告标准实验报告（实验）课程名称学生姓名：李培睿学号：2013020904026指导教师：程建一、实验名称《MATLAB与数值分析》第一次上机实验二、实验目的1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。
（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序）2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。
（用.m 文件编写进行符号因式分解和函数求反的程序）3. 掌握Matlab函数的编写规范。
4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。
（用.m 文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释）5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。
三、实验内容1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。
并以x，y为坐标显示图像x(n+1) = a*x(n)-b*(y(n)-x(n) ); y(n+1) = b*x(n)+a*(y(n)-x(n) )2. 编程实现奥运5环图，允许用户输入环的直径。
3. 实现对输入任意长度向量元素的冒泡排序的升序排列。
不允许使用sort函数。
四、实验数据及结果分析题目一：①在Editor窗口编写函数代码如下：并将编写的函数文件用“draw.m”储存在指定地址；②在Command窗口输入如下命令：③得到图形结果如下：题目二：①在Editor窗口编写函数代码如下：并将编写的函数文件用“circle.m”储存在指定地址；②再次在Editor窗口编写代码：并将编写的函数文件用“Olympic.m”储存在指定地址；③在Command窗口输入如下指令（半径可任意输入）：④按回车执行，将在图形窗口获得五环旗：题目三：①在Editor窗口编写函数代码如下：并用.将编写的函数文件用“qipaofa.m”储存在指定地址；②在Command窗口输入一组乱序数值，则可以得到升序排序结果如下：五、总结及心得体会1. 要熟悉MATLAB编译软件的使用方法，明白有关语法，语句的基本用法，才可以在编写程序的时候游刃有余，不至于寸步难行。
MATLAB中的偏微分方程数值解法
MATLAB中的偏微分方程数值解法偏微分方程（Partial Differential Equations，PDEs）是数学中的重要概念，广泛应用于物理学、工程学、经济学等领域。
解决偏微分方程的精确解往往非常困难，因此数值方法成为求解这类问题的有效途径。
而在MATLAB中，有丰富的数值解法可供选择。
本文将介绍MATLAB中几种常见的偏微分方程数值解法，并通过具体案例加深对其应用的理解。
一、有限差分法（Finite Difference Method）有限差分法是最为经典和常用的偏微分方程数值解法之一。
它将偏微分方程的导数转化为差分方程，通过离散化空间和时间上的变量，将连续问题转化为离散问题。
在MATLAB中，使用有限差分法可以比较容易地实现对偏微分方程的数值求解。
例如，考虑一维热传导方程（Heat Equation）：∂u/∂t = k * ∂²u/∂x²其中，u为温度分布随时间和空间的变化，k为热传导系数。
假设初始条件为一段长度为L的棒子上的温度分布，边界条件可以是固定温度、热交换等。
有限差分法可以将空间离散化为N个节点，时间离散化为M个时刻。
我们可以使用中心差分近似来计算二阶空间导数，从而得到以下差分方程：u(i,j+1) = u(i,j) + Δt * (k * (u(i+1,j) - 2 * u(i,j) + u(i-1,j))/Δx²)其中，i表示空间节点，j表示时间步。
Δt和Δx分别为时间和空间步长。
通过逐步迭代更新节点的温度值，我们可以得到整个时间范围内的温度分布。
而MATLAB提供的矩阵计算功能，可以大大简化有限差分法的实现过程。
二、有限元法（Finite Element Method）有限元法是另一种常用的偏微分方程数值解法，特点是适用于复杂的几何形状和边界条件。
它将求解区域离散化为多个小单元，通过构建并求解代数方程组来逼近连续问题。
在MATLAB中，我们可以使用Partial Differential Equation Toolbox提供的函数进行有限元法求解。
数值分析matlab实验报告
数值分析matlab实验报告数值分析 Matlab 实验报告一、实验目的数值分析是研究各种数学问题数值解法的学科，Matlab 则是一款功能强大的科学计算软件。
本次实验旨在通过使用 Matlab 解决一系列数值分析问题，加深对数值分析方法的理解和应用能力，掌握数值计算中的误差分析、数值逼近、数值积分与数值微分等基本概念和方法，并培养运用计算机解决实际数学问题的能力。
二、实验内容（一）误差分析在数值计算中，误差是不可避免的。
通过对给定函数进行计算，分析截断误差和舍入误差的影响。
例如，计算函数＄f(x) ＝＼sin(x)＄在＄x ＝ 05$ 附近的值，比较不同精度下的结果差异。
（二）数值逼近1、多项式插值使用拉格朗日插值法和牛顿插值法对给定的数据点进行插值，得到拟合多项式，并分析其误差。
2、曲线拟合采用最小二乘法对给定的数据进行线性和非线性曲线拟合，如多项式曲线拟合和指数曲线拟合。
（三）数值积分1、牛顿柯特斯公式实现梯形公式、辛普森公式和柯特斯公式，计算给定函数在特定区间上的积分值，并分析误差。
2、高斯求积公式使用高斯勒让德求积公式计算积分，比较其精度与牛顿柯特斯公式的差异。
（四）数值微分利用差商公式计算函数的数值导数，分析步长对结果的影响，探讨如何选择合适的步长以提高精度。
三、实验步骤（一）误差分析1、定义函数｀compute_sin_error` 来计算不同精度下的正弦函数值和误差。
｀｀｀matlabfunction value, error ＝ compute_sin_error(x, precision)true_value ＝ sin(x)；computed_value ＝ vpa(sin(x)， precision)；error ＝ abs(true_value computed_value)；end｀｀｀2、在主程序中调用该函数，分别设置不同的精度进行计算和分析。
（二）数值逼近1、拉格朗日插值法｀｀｀matlabfunction L ＝ lagrange_interpolation(x, y, xi)n ＝ length(x)；L ＝ 0;for i ＝ 1:nli ＝ 1;for j ＝ 1:nif j ～＝ ili ＝ li （xi x(j)）／（x(i) x(j)）；endendL ＝ L ＋ y(i) li;endend｀｀｀2、牛顿插值法｀｀｀matlabfunction N ＝ newton_interpolation(x, y, xi)n ＝ length(x)；％计算差商表D ＝ zeros(n, n)；D(：， 1) ＝ y'；for j ＝ 2:nfor i ＝ j:nD(i, j) ＝（D(i, j 1) D(i 1, j 1)）／（x(i) x(i j ＋ 1)）；endend％计算插值结果N ＝ D(1, 1)；term ＝ 1;for i ＝ 2:nterm ＝ term （xi x(i 1)）；N ＝ N ＋ D(i, i) term;endend｀｀｀3、曲线拟合｀｀｀matlab％线性最小二乘拟合p ＝ polyfit(x, y, 1)；y_fit_linear ＝ polyval(p, x)；％多项式曲线拟合p ＝ polyfit(x, y, n)；％ n 为多项式的次数y_fit_poly ＝ polyval(p, x)；％指数曲线拟合p ＝ fit(x, y, ＇exp1'）；y_fit_exp ＝ p(x)；｀｀｀（三）数值积分1、梯形公式｀｀｀matlabfunction T ＝ trapezoidal_rule(f, a, b, n)h ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；T ＝ h （（y(1) ＋ y(end)）／ 2 ＋ sum(y(2:end 1)））；end｀｀｀2、辛普森公式｀｀｀matlabfunction S ＝ simpson_rule(f, a, b, n)if mod(n, 2) ～＝ 0error(＇n 必须为偶数'）；endh ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；S ＝ h ／ 3 （y(1) ＋ 4 sum(y(2:2:end 1)）＋ 2 sum(y(3:2:end 2)）＋ y(end)）；end｀｀｀3、柯特斯公式｀｀｀matlabfunction C ＝ cotes_rule(f, a, b, n)h ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；w ＝ 7, 32, 12, 32, 7 ／ 90;C ＝ h sum(w y)；end｀｀｀4、高斯勒让德求积公式｀｀｀matlabfunction G ＝ gauss_legendre_integration(f, a, b)x, w ＝ gauss_legendre(5)；％选择适当的节点数t ＝（b a) ／ 2 x ＋（a ＋ b) ／ 2;G ＝（b a) ／ 2 sum(w f(t)）；end｀｀｀（四）数值微分｀｀｀matlabfunction dydx ＝ numerical_derivative(f, x, h)dydx ＝（f(x ＋ h) f(x h)）／（2 h)；end｀｀｀四、实验结果与分析（一）误差分析通过不同精度的计算，发现随着精度的提高，误差逐渐减小，但计算时间也相应增加。
数值分析在生活中的应用举例及Matlab实现
一、最小二乘法，用MATLAB实现1. 数值实例下面给定的是乌鲁木齐最近1个月早晨7：00左右（新疆时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。
下面用MATLAB编程对上述数据进行最小二乘拟合。
下面用MATLAB编程对上述数据进行最小二乘拟合2、程序代码x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%3、数值结果不同次数多项式拟合误差平方和为：r1=67.6659r2=20.1060r3=3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
4、拟合曲线如下图二、 线性方程组的求解（ 高斯-塞德尔迭代算法 ）1、实例： 求解线性方程组（见书P233页）⎪⎪⎩⎪⎪⎨⎧=++=-+=+-3612363311420238321321321x x x x x x x x x 记A x=b, 其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=363320,,12361114238321b x A x x x任取初始值()()Tx0000=，进行迭代。
复合形法matlab程序及例题
复合形法matlab程序及例题复合形法是一种数值分析方法，用于求解多元函数的最小值。
本文将介绍复合形法的基本原理及其在MATLAB中的实现，同时给出一个具体的例题进行演示。
一、复合形法基本原理复合形法（Nelder-Mead算法）是一种无导数优化方法，它的基本思想是在函数域内随机选择一些点，通过对这些点进行逐步变换、调整，最终找到极小值点。
其算法流程如下：1. 选定n+1个点作为初始点，其中n为被优化函数的自变量个数。
这些点可以随机生成，也可以通过一定的规则生成。
2. 对这些点进行排序，确定最小值点、最大值点和次大值点。
3. 通过一系列变换操作，将最大值点向中间靠拢，同时保持其他点的位置不变。
4. 根据变换后的结果，重新排序，判断是否满足终止条件。
若不满足，继续进行变换操作，直到满足终止条件为止。
终止条件可以是迭代次数、函数值的变化量或某些停止准则的满足等。
二、复合形法MATLAB程序实现复合形法的MATLAB程序可以通过fminsearch函数实现，其语法格式如下：[x,fval,exitflag,output] = fminsearch(fun,x0,options) 其中，fun为被优化的函数，x0为初始点，options为优化选项。
程序返回变量x为最小值点，fval为最小值，exitflag为程序停止状态，output为程序运行信息。
下面是一个简单的例子，用来演示如何使用MATLAB实现复合形法：1. 定义被优化的函数：f(x,y) = (x-2)^2+(y+1)^2+2function z = myfun(x)z = (x(1)-2)^2 + (x(2)+1)^2 + 2;2. 定义初始点和优化选项：x0 = [-1,1];opts = optimset('Display','iter');3. 运行fminsearch函数：[x,fval,exitflag,output] = fminsearch(@myfun,x0,opts);4. 输出结果：x = 1.9999 -0.9999fval = 2.0000exitflag = 1output =iterations: 30funcCount: 58algorithm: 'Nelder-Mead simplex direct search'message: 'Optimization terminated: relative function value decreasing below OPTIONS.TolFun.'最终输出结果表明，通过30次迭代，算法找到了函数的最小值点为x=1.9999，y=-0.9999，函数值为2.0000。
matlab 数值解方程
MATLAB 是一个强大的数值计算和数据分析工具，可以用来求解各种数学方程。
以下是一些基本的 MATLAB 函数，用于求解不同类型的方程。
1. **符号方程求解**：使用 `syms` 和 `solve` 函数。
```matlab
syms x
solution = solve(x^2 + 3*x - 4 = 0);
```
2. **数值方程求解**：使用 `fzero` 函数。
```matlab
f = @(x) x^2 - 3; % 定义函数
solution = fzero(f, 1); % 在区间 [1,2] 内寻找零点
```
3. **系统方程求解**：使用 `fsolve` 函数。
```matlab
fun = @(x) [x(1)^2 + x(2)^2 - 4; x(1) + x(2) - 2]; % 定义系统函数
solution = fsolve(fun, [1,1]); % 使用初始值 [1,1] 进行求解
```
4. **非线性方程求解**：使用 `fminbnd` 或 `fminunc` 函数。
```matlab
fun = @(x) x^3 - x - 1; % 定义函数
solution = fminbnd(fun, -10, 10); % 在区间 [-10,10] 内寻找最小值点
```
注意：以上所有的解决方案都需要先定义方程或系统，然后选择合适的初始值和参数范围进行求解。
MATLAB 的文档提供了更详细的例子和说明，可以帮助你理解如何使用这些函数。
分数阶微分方程的数值解法及其MATLAB实现
分数阶微分方程的数值解法及其MATLAB实现1.分数阶微分方程的定义和性质分数阶导数是一般导数的推广，可以用分数阶微分方程来描述分数阶导数的性质。
分数阶微分方程一般形式如下：D^αy(t)=f(t,y(t)),0<α≤1其中，α是分数阶指数，D^α是分数阶导数符号，y(t)是未知函数，f(t,y(t))是已知的函数。
2.分数阶微分方程的数值解法由于分数阶导数的定义较为复杂，常见的数值解法有两种：格点差分法和Laplace变换法。
2.1格点差分法格点差分法是将连续的分数阶导数问题离散化为离散的整数阶微分方程问题，然后利用传统的整数阶微分方程数值解法求解。
具体步骤如下：(1)将时间区间[0,T]平均分为N段，使得Δt=T/N。
(2)将分数阶导数D^α近似为整数阶导数D_m^β，其中m>α>m-1，β=m-α。
(3)将方程D^αy(t)=f(t,y(t))离散为差分方程D_m^βy(t_k)≈f(t_k,y(t_k))，其中t_k=kΔt，k=0,1,2,...,N。
(4)解差分方程，得到数值解y_k，k=0,1,2,...,N。
2.2 Laplace变换法Laplace变换法是将分数阶微分方程转化为Laplace变换形式的整数阶微分方程，然后利用Laplace变换的性质和经典的整数阶微分方程数值解法求解。
具体步骤如下：(1) 对分数阶微分方程D^αy(t) = f(t, y(t)) 进行Laplace变换，得到整数阶微分方程s^αY(s) - s^αy(0) + s^α-1Y(s) = F(s)，其中Y(s) 和 F(s) 分别为 y(t) 和 f(t,y(t)) 的Laplace变换。
(2)将Y(s)和F(s)用有限差分近似替代，并将方程离散化为差分方程，得到Y_k(s)和F_k(s)，其中k是离散时间步长。
(3) 解差分方程，得到 Y_k。
将 Y_k 用逆Laplace变换转换为 y_k。
MATLAB 数值分析
13.2
极小化
作图除了提供视觉信息外，还常常需要确定一个函数的其它更多的特殊属性。在许多 应用中，特别感兴趣的是确定函数的极值，即最大值（峰值）和最小值（谷值）。数学上， 可通过确定函数导数（斜率）为零的点，解析上求出这些极值点。检验 humps 的图形在峰 值和谷值点上的斜率就很容易理解这个事实。显然，如果定义的函数简单，则这种方法常 常奏效。然而，即使很多容易求导的函数，也常常很难找到导数为零的点。在这种情况下， 以及很难或不可能解析上求得导数的情况下，必须数值上寻找函数的极值点。MATLAB 提 供了两个完成此功能的函数 fmin 和 fmins。 这两个函数分别寻找一维或 n 维函数的最小值。 这里仅讨论 fmin。有关 fmins 的详细信息，参阅《MATLAB 参考指南》。因为 f(x)的最大 值等于-f(x)的最小值，所以，上述 fmin 和 fmins 可用来求最大值和最小值。如果还不清楚， 把上述图形倒过来看，在这个状态下，峰值变成了谷值，而谷值则变成了峰值。 为了解释求解一维函数的最小值和最大值， 再考虑上述例子。 从图 13.2 可知， 在 xmax=0.7 附 近 有 一 个 最 大 值 ， 并 且 在 xmin=4 附 近 有 一 个 最 小 值 。 而 这 些 点 的 解 析 值 为 ： x m a x / 4 0.785 和 x min 5 / 4 393 . 。为了方便，用文本编辑器编写一个脚本 M 文件，并用 fmin 寻出数值上极值点，给出函数主体如下： % ex_fmin.m fn=‘ 2*exp(-x)*sin(x) ‘; xmin=fmin(fn , 2 , 5)
% better approximation
自然地，上述两个结果不同。基于对图形的观察，粗略近似可能低估了实际面积。除 非特别精确，没有准则说明哪种近似效果更好。很明显，如果人们能够以某种方式改变单 个梯形的宽度，以适应函数的特性，即当函数变化快时，使得梯形的宽度变窄，这样就能 够得到更精确的结果。 MATLAB 的函数 quad 和 quad8 是基于数学上的正方形概念来计算函数的面积， 这些 积分函数的操作方式一样。为获得更准确的结果，两个函数在所需的区间都要计算被积函 数。 此外， 与简单的梯形比较， 这两个函数进行更高阶的近似， 而且 quad8 比 quad 更精确。 这两个函数的调用方法与 fzero 相同，即 >>area=quad(‘ humps ‘ , -1 , 2) % find area between -1 and 2
数值分析matlab方法插值法
其中,
n
【注】
x [a, b] w(x) (x x j ) j 0
（1）误差估计
Rn (x)
M n1 (n 1)!
w( x)
M n1
max
x( a ,b )
f
(n1) (x)
（2）余项与 x、M n1 节点的位置、个数 n 有关
（3）当 f (x)是 n 的多项式时Ln (x) f (x) n
M2 2!
(x
x0 )(x
x1 )
其中，
M2
max
x( x0 , x1 )
f
(x)
x
[
6
,
4
]
，所以
R1
(
5
24
)
sin
2!
4 (5
24
)( 5
6 24
)
4
0.0061
2）
抛物插值误差估计.因为
R2 (x)
M3 3!
(x
x0 )(x
x1)(x
x2 )
其中，
M3
max
x( x0 ,x2 )
f (x)
yiynewtonbackwardxyxicos035yi09394数值分析插值法55埃尔米特插值2n12n2数值分析插值法551埃尔米特插值多项式的存在唯一性2n1数值分析插值法数值分析插值法552埃尔米特插值余项553三次埃尔米特插值多项式maxsinxsin1数值分析插值法56561高次插值的病态性质0908原函数150706050405分段线性插值0302543214321055数值分析插值法562分段低次插值方法563分段低次插值余项090807060504030201090807060504030201543214321数值分析插值法57三次样条插值571三次样条插值572三弯矩法数值分析插值法573三次样条插值的误差估计与收敛性58插值运算的matlab函数581一维插值函数interp1yiinterp1xyximethod?linear?yiinterp1xyxilinear?1200135019
matlab算法-求解微分方程数值解和解析解
MATLAB是一种用于数学计算、工程和科学应用程序开发的高级技术计算语言和交互式环境。
它被广泛应用于各种领域，尤其在工程和科学领域中被用于解决复杂的数学问题。
微分方程是许多工程和科学问题的基本数学描述，求解微分方程的数值解和解析解是MATLAB算法的一个重要应用。
1. 求解微分方程数值解在MATLAB中，可以使用各种数值方法来求解微分方程的数值解。
其中，常见的方法包括欧拉法、改进的欧拉法、四阶龙格-库塔法等。
这些数值方法可以通过编写MATLAB脚本来实现，从而得到微分方程的近似数值解。
以常微分方程为例，可以使用ode45函数来求解微分方程的数值解。
该函数是MATLAB中用于求解常微分方程初值问题的快速、鲁棒的数值方法，可以有效地得到微分方程的数值解。
2. 求解微分方程解析解除了求解微分方程的数值解外，MATLAB还可以用于求解微分方程的解析解。
对于一些特定类型的微分方程，可以使用符号计算工具箱中的函数来求解微分方程的解析解。
通过符号计算工具箱，可以对微分方程进行符号化处理，从而得到微分方程的解析解。
这对于研究微分方程的性质和特点非常有帮助，也有助于理论分析和验证数值解的准确性。
3. MATLAB算法应用举例在实际工程和科学应用中，MATLAB算法求解微分方程问题非常常见。
在控制系统设计中，经常需要对系统的动态特性进行分析和设计，这通常涉及到微分方程的建模和求解。
通过MATLAB算法，可以对系统的微分方程进行数值求解，从而得到系统的响应曲线和动态特性。
另外，在物理学、生物学、经济学等领域的建模和仿真中，也经常需要用到MATLAB算法来求解微分方程问题。
4. MATLAB算法优势相比于其他数学软件和编程语言，MATLAB在求解微分方程问题上具有明显的优势。
MATLAB提供了丰富的数值方法和工具，能够方便地对各种微分方程进行数值求解。
MATLAB具有直观的交互式界面和强大的绘图功能，能够直观地展示微分方程的数值解和解析解，有利于分析和理解问题。
数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。
由角动量定理我们知道εJ M =化简得到0sin 22=+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。
实际上这是一个解二阶常微分方程的问题。
θ在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程W T =h mg ∆=2J 21ω化简得到四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。
所以比欧拉稳定。
运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。
通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。
三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0; %初值%RK_z(1)=0;初值for i=1:length(x)-1K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1 and L1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;% K4L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variable x');ylabel('Variable y');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')%周期fprintf('周期等于%d',T*h)两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;%欧拉的初值Euler_z(1)=0;RK_y(1)=0;%龙格库塔初值RK_z(1)=0;for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 andL1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);% K2 and L2K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);% K3 and L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);%角度的峰值也就是πfprintf('角度峰值等于%d',y)fprintf('\n')%周期fprintf('周期等于%d',T*h)xlabel('时间');ylabel('角度');legend('欧拉','RK4');。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第2章牛顿插值法实现
参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55.
求牛顿插值多项式和差商的MA TLAB 主程序：
function[A,C,L,wcgs,Cw]=newpoly(X,Y)
n=length(X);A=zeros(n,n);A(:,1) =Y';
s=0.0;p=1.0;q=1.0;c1=1.0;
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));
end
b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1;
end
C=A(n,n);b=poly(X(n));q1=conv(q1,b);
for k=(n-1):-1:1
C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k);
end
L(k,:)=poly2sym(C);Q=poly2sym(q1);
syms M
wcgs=M*Q/c1;Cw=q1/c1;
（1）保存名为newpoly.m 的M 文件
（2）输入MA TLAB 程序
>> X=[242,243,249,250];
>> Y=[13.681,13.526,13.098,13.095];
>> [A,C,L,wcgs,Cw]=newpoly(X,Y)
输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其
)
()()1(ξ+n n f x R 的系数向量Cw 。
A =
13.6810 0 0 0
13.5260 -0.1550 0 0
13.0980 -0.0713 0.0120 0
13.0950 -0.0030 0.0098 -0.0003
C =
1.0e+003 *
-0.0000 0.0002 -0.0551 4.7634
L =
- (23*x^3)/84000 + (2981*x^2)/14000 - (7757472138947345*x)/140737488355328 + 5237382665812919/1099511627776
wcgs =
(M*(x^4 - 984*x^3 + 363071*x^2 - 59535444*x + 3660673500))/24
Cw =
1.0e+008 *
0.0000 -0.0000 0.0002 -0.0248 1.5253
输入MATLAB程序
>> x=244;
>> y=- (23*x^3)/84000 + (2981*x^2)/14000 - (7757472138947345*x)/140737488355328 + 5237382665812919/1099511627776
y =
13.3976
输入MATLAB程序
>> x=[244,245,246,247,248];
>> y=- (23.*x.^3)./84000 + (2981.*x.^2)./14000 - (7757472138947345.*x)./140737488355328 + 5237382665812919./1099511627776
y =
13.3976 13.2943 13.2143 13.1560 13.1178
第4章 高斯-勒让德积分公式实现
用高斯-勒让德积分公式计算dx e I x ⎰--=112221
π，取代数精度为3和5，再根据截断误差公
式写出误差公式，并将计算结果与精确值进行比较。
（1）建立并保存fun.m 文件命名的M 文件函数
function y=fun(x)
y=exp((-x.^2)./2)./(sqrt(2*pi));
（2）建立并保存GaussR1.m 文件命名的M 文件函数
function[GL,Y,RGn]=GaussR1(fun,X,A)
n=length(X);n2=2*n;Y=feval(fun,X);GL=sum(A.*Y);
sun=1;su2n=1;su2n1=1;wome=1;
syms x
for k=1:n
wome=wome*(x-X(k));
end
wome2=wome^2;Fr=int(wome2,x,-1,1);
for k=1:n2
su2n=su2n*k;
end
syms M
RGn=Fr*M/su2n;
（3）输入程序
代数精度为3：
>> X1=[-1/(3^(1/2)),1/(3^(1/2))];A1=[1,1];
>> [GL1,Y1,Rn1]=GaussR1(@fun,X1,A1)
GL1 =
0.6754
Y1 =
0.3377 0.3377
Rn1 =
M/135
代数精度为5：
>> X2=[-(3/5)^(1/2),0,(3/5)^(1/2)];A2=[5/9,8/9,5/9];
>> [GL2,Y1,Rn2]=GaussR1(@fun,X2,A2)
GL2 =
0.6830
Y1 =
0.2955 0.3989 0.2955
Rn2 =
M/15750
截断误差：
>> syms x
>> fi=int(exp((-x.^2)./2)./(sqrt(2*pi)),x,-1,1);
>> Fs=double(fi),wGL1=double(abs(fi-GL1)),wGL2=double(abs(fi-GL2)) Fs =
0.6827
wGL1 =
0.0073
wGL2 =
3.0777e-004。