工程数值分析实验(龙格库塔,最小二乘法)
计算方法上机作业——龙格库塔法
计算方法上机作业——龙格库塔法龙格库塔法（Runge-Kutta method）是一种常用于求解常微分方程（Ordinary Differential Equation，ODE）的数值解法。
它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）在20世纪初提出的。
龙格库塔法的基本思想是通过数值逼近来计算微分方程的近似解。
在讲解龙格库塔法之前，我们先来简单回顾一下ODE的一阶常微分方程的基本形式：y′(y)=y(y,y)其中，y(y,y)是已知函数。
龙格库塔法的核心是使用差分逼近计算函数的斜率。
假设我们要求解的方程为：y′(y)=y(y,y),y(y)=y₀所需计算的点为y₀,y₁,...,yy，对应的函数值为y₀,y₁,...,yy，其中y是步长的个数。
龙格库塔法通过递推关系式来计算估计值，并不断更新当前点的函数值。
接下来以龙格库塔法的经典四阶形式为例进行说明。
该方法的基本方程如下：yy+1=yy+(y₁+2y₂+2y₃+y₄)/6y₁=ℎy(yy,yy)y₂=ℎy(yy+ℎ/2,yy+y₁/2)y₃=ℎy(yy+ℎ/2,yy+y₂/2)y₄=ℎy(yy+ℎ,yy+y₃)其中y表示当前步骤，ℎ表示步长，yy表示当前点的函数值，y₁,y₂,y₃和y₄则表示对应的斜率。
使用龙格库塔法，我们可以通过不断递归计算来求得指定区间（例如[y,y]）上的函数值。
具体步骤如下：1.确定求解区间[y,y]和初始点(y₀,y₀)以及步长ℎ。
2.初始化：设置yy=y₀，yy=y₀。
3.对所有y=0,...,y−1：计算y₁,y₂,y₃和y₄，根据上述递推关系式。
根据递推关系式计算yy+1更新当前点的函数值，即yy+1=y(yy+1)。
更新当前点的y值，即yy+1=yy+ℎ。
4.返回结果：最终求得的函数值。
需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。
龙格库塔
数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
经典四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。
令初值问题表述如下。
则，对于该问题的RK4由如下方程给出：其中这样，下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。
该斜率是以下斜率的加权平均：∙k1是时间段开始时的斜率；∙k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点t n + h/2的值；∙k3也是中点的斜率，但是这次采用斜率k2决定y值；∙k4是时间段终点的斜率，其y值用k3决定。
当四个斜率取平均时，中点的斜率有更大的权值：RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。
注意上述公式对于标量或者向量函数(y可以是向量)都适用。
显式龙格库塔法显示龙格－库塔法是上述RK4法的一个推广。
它由下式给出其中如果要求方法有精度p则还有相应的条件，也就是要求舍入误差为O(h p+1)时的条件。
这些可以从舍入误差本身的定义中导出。
例如，一个2阶精度的2段方法要求b1 + b2 = 1, b2c2 = 1/2, 以及b2a21 = 1/2。
在Matlab下输入：edit，然后将下面两行百分号之间的内容，复制进去，保存%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%function dxdt=ode_Miss_ghost(t,x)%分别用x(1),x(2),x(3),x(4)代替N1,P1,N2,P2N1=x(1);P1=x(2);N2=x(3);P2=x(4);K=2;tau_c=3e-9;tan_p=6e-12;beta =5e-5;delta=0.692;eta =0.0001;fm =8e6;Ith =26e-3;Ib =1.5*Ith;Im =0.3*Ith;I1=Ib+Im*sin(2*pi*fm*t)+K*P2;I2=Ib+Im*sin(2*pi*fm*t)+K*P1;dxdt=[(I1/Ith-N1-(N1-delta)/(1-delta)*P1)/tau_e;((N1-delta)/(1-delta)*(1-eta*P1)*P1-P1+beta*N1)/tau_p;(I2/Ith-N2-(N2-delta)/(1-delta)*P2)/tau_e;((N2-delta)/(1-delta)*(1-eta*P2)*P2-P2+beta*N2)/tau_p;]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%在Matlab下面输入：t_start=0;t_end=2e-9;y0=[1e-3;1e-4;0;0]; %初值[x,y]=ode15s('ode_Miss_ghost',[0,t_end],y0);plot(x,y);legend('N1','P1','N2','P2');xlabel('x');Mathlab定步长龙格－库塔法[345阶]、自适应步长rkf45经常看到很多朋友问定步长的龙格库塔法设置问题，下面吧定步长三阶、四阶、五阶龙格库塔程序贴出来，有需要的可以看看ODE3 三阶龙格－库塔法CODE:function Y = ode3(odefun,tspan,y0,varargin)%ODE3 Solve differential equations with a non-adaptive method of order 3.% Y = ODE3(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates% the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE3(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN% but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the Bogacki-Shampine Runge-Kutta method of order 3.%% Example% tspan = 0:0.1:20;% y = ode3(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1, % and plots the first component of the solution.%if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);F = zeros(neq,3);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);F(:,1) = feval(odefun,ti,yi,varargin{:});F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});F(:,3) = feval(odefun,ti+0.75*hi,yi+0.75*hi*F(:,2),varargin{:});Y(:,i) = yi + (hi/9)*(2*F(:,1) + 3*F(:,2) + 4*F(:,3));endY = Y.';ODE4 四阶龙格－库塔法CODE:function Y = ode4(odefun,tspan,y0,varargin)%ODE4 Solve differential equations with a non-adaptive method of order 4.% Y = ODE4(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE4(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the classical Runge-Kutta method of order 4.%% Example% tspan = 0:0.1:20;% y = ode4(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1,% and plots the first component of the solution.%if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr]; error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);F = zeros(neq,4);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);F(:,1) = feval(odefun,ti,yi,varargin{:});F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:}); F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:}); F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));endY = Y.';定步长RK4(自编):CODE:function Y=RungeKutta4(f,xn,y0)% xn=0:.1:1;% y0=1;% y_n=[];% f=@(X1,Y1) Y1-2*X1/Y1;y_n=[];h=diff(xn(1:2));for i=1:length(xn)-1K1=f(xn(i),y0);K2=f(xn(i)+h/2,y0+h*K1/2);K3=f(xn(i)+h/2,y0+h*K2/2);K4=f(xn(i)+h,y0+h*K3);y_n=[y_n;y0+h/6*(K1+2*K2+2*K3+K4)];y0=y_n(end);endY=y_n;ODE5 五阶龙格－库塔法CODE:function Y = ode5(odefun,tspan,y0,varargin)%ODE5 Solve differential equations with a non-adaptive method of order 5.% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the Dormand-Prince method of order 5 in a general % framework of explicit Runge-Kutta methods.%% Example% tspan = 0:0.1:20;% y = ode5(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1,% and plots the first component of the solution.if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);% Method coefficients -- Butcher's tableau%% C | A% --+---% | BC = [1/5; 3/10; 4/5; 8/9; 1];A = [ 1/5, 0, 0, 0, 03/40, 9/40, 0, 0, 044/45 -56/15, 32/9, 0, 019372/6561, -25360/2187, 64448/6561, -212/729, 0 9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]; % More convenient storageA = A.';B = B(:);nstages = length(B);F = zeros(neq,nstages);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);% General explicit Runge-Kutta frameworkF(:,1) = feval(odefun,ti,yi,varargin{:});for stage = 2:nstageststage = ti + C(stage-1)*hi;ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));F(:,stage) = feval(odefun,tstage,ystage,varargin{:});endY(:,i) = yi + F*(hi*B);endY = Y.';-------------------------------------------------------------------------------------------------------------------自适应步长RKF45(相当于ode45)ODE45 是4阶方法提供候选解，5阶方法控制误差。
数值分析实验报告
数值分析实验报告
一、实验背景
本实验主要介绍了数值分析的各种方法。
在科学计算中，为了求解一
组常微分方程或一些极限问题，数值分析是一种有用的方法。
数值分析是
一种运用计算机技术对复杂模型的问题进行数学分析的重要手段，它利用
数学模型和计算机程序来解决复杂的数学和科学问题。
二、实验内容
本实验通过MATLAB软件，展示了以下几种数值分析方法：
（1）拉格朗日插值法：拉格朗日插值法是由法国数学家拉格朗日发
明的一种插值方法，它可以用来插值一组数据，我们使用拉格朗日插值法
对给定的点进行插值，得到相应的拉格朗日多项式，从而计算出任意一个
点的函数值。
（2）最小二乘法：最小二乘法是一种常用的数据拟合方法，它可以
用来拟合满足一定函数的点的数据，它的主要思想是使得数据点到拟合曲
线之间的距离的平方和最小。
（3）牛顿插值法：牛顿插值法是一种基于差商的插值方法，它可以
用来插值一组数据，可以求得一组数据的插值函数。
（4）三次样条插值：三次样条插值是一种基于三次样条的插值方法，它可以用来对一组数据进行插值，可以求得一组数据的插值函数。
三、实验步骤
1．首先启动MATLAB软件。
四阶龙格库塔实验报告
三、四阶Runge-Kutta 法求解常微分方程一、龙格库塔法的思想根据第九章的知识可知道，Euler 方法的局部截断误差是2()O h ，而当用Euler 方法估计出1,()(1)n n n n y y hf x y +=+ 再用梯形公式111[(,)(,)](2)2n n n n n n h y y f x y f x y +++=++进行校正，即采用改进Euler 方法得出数值解的截断误差为3()O h 。
由Lagrange 微分中值定理'11()()()()()(,())(3)n n n n n y x y x y x x y x hf y ξξξ++=+-=+ 记*(,())k hf y ξξ=，得到*1()()(4)n n y x y x k +=+这样只要给出一种计算*k 的算法，就能得到相应的计算公式。
用这种观点的来分析Euler 方法和改进Euler 方法，Euler 方法的迭代公式可改写为111(,)n n n n y y k k hf x y +=+=改进Euler 方法的预报-校正公式可改写为 1121211()2(,),(,)n n n n n n y y k k k hf x y k hf x h y k +=++==++ Euler 方法实际上是用一个点处的值1k 近似*k ，而改进Euler 方法是用两个点处的值1k ，和2k ，做算术平均值近似*k 自然改进Euler 方法要优于Euler 方法。
因此，可以想到假如在1[,]n n x x +内多预报几个点值i k ，并用他们的加权平均值作为*k 的近似值，则有可能构造出具有更高精度的计算公式，这就是Runge-Kutta 法的基本思想。
二、四阶龙格库塔法由Runge-Kutta 的基本思想，构造四阶Runge-Kutta 法是利用1234,,k k k k 和的加权平均值来近似*k ，因此令1112233441211132211243312213(,)(,)(5)(,)(,)n n n n n n n n n n y y w K w K w K w K K hf x y K hf x h y K K hf x h y K K K hf x h y K K K αβαβγαβγη+=++++⎧⎪=⎪⎪=++⎨⎪=+++⎪⎪=++++⎩使得511()()n n y x y O h ++-=即其总体截断误差为4()O h 。
数值计算方法第3、4次实验--龙贝格--龙格库塔
数值计算方法第3、4次实验--龙贝格--龙格库塔完成复合梯形、复合辛普森求积公式的编写，使用变步长和事后误差估计的方法计算P145例题6.3.1对应表6-3，课后作业题7.以下是变步长梯形求积法的MATLAB实现：___(a,b,eps)m=1;t0=0;t2=(f(a)+f(b))/4+f((a+b)/2);while(abs(t2-t0)>eps)m=m+1;p=t2;t1=0;n=2^m;h=(b-a)/n;for(k=0:(n-1))t1=t1+h*((f(a+k*h)+f(a+(k+1)*h))/4+f(a+(k+1/2)*h)/2); endt2=t1;t0=p;endI=t2;接下来是对该函数的调用示例：AdaptiveTrapezoidal(1,2,0.)m =1t2 =3.xxxxxxxxm =2t2 =2.xxxxxxxxm =3t2 =2.xxxxxxxxm =4t2 =2.xxxxxxxxm =5t2 =2.xxxxxxxxm =6t2 =2.xxxxxxxxm =7t2 =2.xxxxxxxxm =8t2 =2.xxxxxxxxm =9t2 =2.xxxxxxxxI =2.xxxxxxxx再看下面的调用示例：AdaptiveTrapezoidal(1,3,0.0001) m =1t2 =7.xxxxxxxxm =2t2 =10.xxxxxxxxm =3t2 =10.xxxxxxxxm =4t2 =10.xxxxxxxxm =5t2 =10.xxxxxxxxm =6t2 =10.xxxxxxxxm =7t2 =10.xxxxxxxxm =8t2 =10.xxxxxxxxI =10.xxxxxxxx接下来是复合梯形求积公式的MATLAB实现：复合梯形求积公式的MATLAB实现姓名：___学号：xxxxxxxx48经过上述修改和改写后，文章已经没有格式错误和明显的问题段落了。
数值分析曲线拟合的最小二乘法实验报告
数值分析曲线拟合的最小二乘法实验报告数值分析曲线拟合的最小二乘法实验报告篇一：数值分析设计曲线拟合的最小二乘法曲线拟合的最小二乘法一、目的和意义在科学实验的统计方法研究中，往往要从一组实验数据?xi,yi??i?0,1,2,?,m?中，寻找自变量x与因变量y之间的函数关系y?F?x?。
由于观测数据往往不准确，因此不要求y?F?x?经过所有点?xi,yi?，而只要求在给定点xi上误差而只要求所在所有给定点xi上的误差?i?F(xi)?yi ?i?0,1,2,?,m?按某种标准最小。
若记????0,?1,?2,?,?m?，就是要求向量?的范数如果用最大范数，计算上困难较大，通常采用欧式范数?最小。
2T 作为误差度量的标准。
F?x?的函数类型往往与实验的物理背景以及数据的实际分布有关，它一般含有某些待定参数。
如果F?x?是所有待定参数的线性函数，那么相应的问题称为线性最小二乘问题，否则称为非线性最小二乘问题。
最小二乘法还是实验数据参数估计的重要工具。
这是因为这种方法比其他方法更容易理解，即使在其他方法失效的情况下，用最小二乘法还能提供解答，而且从统计学的观点分析，用该方法求得各项估计具有最优统计特征，因此这一方法也是系统识别的重要基础。
线性最小二乘问题可以借助多元微分学知识通过求解法方程组得到解答。
用最小二乘法求拟合曲线时，首先要确定S?x?的形式。
这不单纯是数学问题，还与所研究问题的运动规律以及所得观测数据?xi,yi?有关；通常要从问题的运动规律以及给定数据描图，确定S?x?的形式，并通过实际计算选出较好的结果。
为了使问题的提法更有一般性，通常把最小二乘法中的? 22 都考虑为加权平方和22 ? ????xi???S?xi??f?xi??? i?0 m 2 这里??xi??0是?a,b?上的加权函数，它表示不同点?xi,f?xi?处的数据比重不同。
?二、计算方法在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y与时间t的拟合曲线。
数值分析实验指导书
《数值分析》实验指导书潍坊学院数学与信息科学学院2012年04月目录目录.............................................................................................. 错误！未定义书签。
实验一插值与曲线拟合的最小二乘法.................................. 错误！未定义书签。
实验二数值积分...................................................................... 错误！未定义书签。
实验三解线性方程组的直接法.............................................. 错误！未定义书签。
实验四解线性方程组的迭代法.............................................. 错误！未定义书签。
实验五非线性方程的数值解法.............................................. 错误！未定义书签。
实验六常微分方程数值解法.................................................... 错误！未定义书签。
实验一 插值与曲线拟合的最小二乘法一、实验目的：1．了解拉格朗日插值法、牛顿插值法、曲线拟合最小二乘法的基本原理和方法；2．掌握拉格朗日插值多项式牛顿插值多项式的用法；3．掌握最小二乘原理，会求拟合函数及超定方程组的最小二乘解。
二、实验内容：1．用拉格朗日插值公式和牛顿插值公式确定函数值； 2．对函数f (x )进行拉格朗日插值和牛顿插值； 3．利用Polyfit 拟合幂函数，利用Polyfit 拟合多项式。
三、实验过程：1．给定函数四个点的数据如下:117.2)1.5(,651.4)9.3(,276.4)3.2(,887.3).11(====f f f f ，试用插值公式确定函数在234.4,101.2=x 处的函数值)(x f 。
第四讲龙格-库塔方法
龙格-库塔方法3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的一般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值方法，若用显式单步法(3.2.2)当，即数值求积用左矩形公式，它就是Euler法(3.1.2)，方法只有一阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶精度的.若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将(3.2.1)右端积分表示为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)一样，用前面已算得的f值表示为(3.2.3)，一般情况可将(3.2.2)的 表示为(3.2.4)其中这里均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K方法.它每步计算r个f值(即)，而k由前面(i-1)个已算出的表示，故公式是显式的.例i如当r=2时，公式可表示为(3.2.5) 其中.改进Euler 法(3.1.11)就是一个二级显式R-K 方法.参数取不同的值，可得到不同公式.3.2.2 二、三级显式R-K 方法对r=2的显式R-K 方法(3.2.5)，要求选择参数，使公式的精度阶p 尽量高，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6) 令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯一.因r=2，由（3.2.5）式可知，于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多，只要都可得到二阶精度R-K方法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)及中点公式(3.2.9)是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求4个待定参数，验证得出是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对r=3的情形，要计算三个k值，即其中将按二元函数在处按Taylor公式展开，然后代入局部截断误差表达式，可得可得三阶方法，其系数共有8个，所应满足的方程为(3.2.10)这是8个未知数6个方程的方程组，解也是不唯一的，通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法：(3.2.11)附：R-K 的三级Kutta 方法程序如下function y = DELGKT3_kuta(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+1K1 = 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 y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K 方法及步长的自动选择利用二元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K 方法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，而33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。
实验龙格—库塔方法
实验 龙格—库塔方法实验目的： 深入理解龙格—库塔方法的原理，同时学会用迭代方法看待某些特定问题。
比较本方法与其它方法（尤其是欧拉方法）解题的不同处，分析两者的精度。
实验要求: 编写程序，用龙格—库塔方法求解第三章3题（修改）实验内容：题目： 取h =0.1，用四阶龙格-库塔方法求解初值问题： 41,211)(22≤≤-+='x y xx f 0)0(=y 并与精确解21x x y +=比较结果. 原理： 本题很明显为迭代方法的一种。
已知初值和函数的导数。
再根据四阶龙格—库塔公式：k 1=dy(x0,y0);k 2=dy((2 * x0 + h)/2 , y0 + h/2 * k 1);k 3=dy((2 * x0 + h)/2 , y0 + h/2 * k 2);k 4=dy(x0 + h , y0 + h * k 3);y0=y0 + h/6 * (k 1 + 2*k 2 + 2*k 3 + k 4);x0=x0+h;设计思想：由四阶龙格—库塔公式很容易看出，只要作一个循环即可完与此项操作。
每次更新k 1,k 2,k 3,k 4,y0以及x0的值。
这样即可得出y 。
源代码：disp('第三章3题：取h=0.1，用四阶龙格-库塔方法求解初值问题：')disp(' diff(y,x,1)=1/(1+x.^2)-2*y.^2 1<=x<=4') disp(' y(0)=0 ')disp(' 并与精确解y=x/(1+x.^2)比较结果.')x0=1; %初值为1y0=0.5; %初值为0.5h=input('请输入步长>>')for i=1:1:5k1=dy(x0,y0);k2=dy((2 * x0 + h)/2 , y0 + h/2 * k1);k3=dy((2 * x0 + h)/2 , y0 + h/2 * k2);k4=dy(x0 + h , y0 + h * k3);y0=y0 + h/6 * (k1 + 2*k2 + 2*k3 + k4);x0=x0+h;fprintf('\nk%d = %f ,k%d = %f ,k%d = %f ,k%d = %f ',1,k1,2,k2,3,k3,4,k4) fprintf('\n =>>y[%d]=%f ',i,y0);endfprintf('\n 以下为用公式：y=x/(1+x.^2)求得的精确值\n')x0=1;for i=0:1:5y0=f(x0);x0=x0+h;fprintf('y[%d]=%f ',i,y0);endfprintf('\n/****************************************************/')(下图为yy=dy(x)以及y=f(x)两函数的源代码：)文件名： f.mfunction y=f(x)y=x./(1+x.^2);文件名：dy.mfunction yy=dy(x,y)yy=1/(1+x.^2)-2*y.^2;实验结果：图形实验体会：本实验中，关于龙格—库塔方法的使用，使我个人又对迭代方法有了新的一种理念。
龙格-库塔方法-文档资料
c3a 2b32
c3a3 1
6
; 2
O (h4)
常见的2种三阶方法：
库塔三阶方法
h yn1yn6(k14k2k3)
k1
f(xn,yn);
k2
hh f(xn2,yn2k1)
k 3 f(x n h ,y n h k 1 2 h k 2 ) •5
四级方法：N = 4
局部截断误差 O ( h 5 )
可见误差随着 x n 的增加呈指数函数递减
当 f y 0 时，微分方程是不稳定的； 而 f y 0 时，微分方程是稳定的。
上面讨论的稳定性，与数值方法和方程中 f 有关
•21
实验方程： y y C ,R e () 0
D e f 3 对单步法 yn 1ynh(xn,yn,h )应用实验方程，
e n 1 e n h [ ( x n ,y ( x n ) , h ) ( x n ,y n , h ) ] T n 1
•15
因为单步法是 p 阶的：h0,0hh0满足|Tn1|Chp1
|e n 1| |e n| h L |e n| C h p 1|en |
其中 1hL,C hp1
•18
三、绝对稳定性 /*Absolute Stibility*/ 计算过程中产生的舍入误差对计算结果的影响
首先以Euler公式为例，来讨论一下舍入误差的传播:
yn1ynhf(xn,yn)
设实际计算得到的点 x n 的近似函数值为 yn yn n，
其中 y n 为精确值， n 为误差
yn1ynhf(xn,yn)
通过适当选取参数 1,2和 p 的值，使得公式具有 2阶精度！！
•3
由泰勒公式展开，要使公式具有 2 阶精度，只需
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
工程数值分析实验报告
指导老师
班级
学号
姓名
实验一：最小二乘法拟合曲线实验
一、实验名称：最小二乘法拟合曲线实验
实验时间： 2015-5-14 实验地点： 主楼机房 实验器材： 计算机matlab
二、实验目的：学会用最小二乘法求拟合数据的多项式，并应用算法于实际问题。
三、实验要求：
（1）根据最小二乘法和加权最小二乘法的基本理论，编写程序构造拟合曲线的法方程，要求可以方便的调整拟合多项式的次数；
（2）采用列主元法解（1）中构造的法方程，给出所拟合的多项式表达式； （3）编写程序计算所拟合多项式的均方误差，并作出离散函数 和拟合函数的图形； （4） 用MATLAB 的内部函数polyfit 求解上面最小二乘法曲线拟合多项式的系数及平方误差，并用MATLAB 的内部函数plot 作出其图形，并与（1）的结果进行比较。
四、算法描述（实验原理与基础理论）
基本原理：从整体上考虑近似函数 同所给数据点 (i=0,1,…,m)误差 (i=0,1,…,m)的大小，常用的方法有以下三种：一是误差 (i=0,1,…,m)绝
对值的最大值
，即误差 向量
的∞—范数；二是误差绝对值的和
，即误差向量r 的1—范数；三是误差平方和 的算术平方根，即误差向量r 的
2—范数；前两种方法简单、自然，但不便于微分运算 ，后一种方法相当于考虑 2—范数的平方，因此在曲线拟合中常采用误差平方和 来 度量误差 (i=0，1，…，m)的整体
大小。
五、实验内容：共有两组给定数据，把给定的数据拟合成多项式。
第一组给定数据点如表1所示如下：
表1 数据表
i
x 0 0．5 0．6 0．7 0．8 0．9 1．0 i
y
1
1．75
1．96
2．19
2．44
2．71
3．00
表2 数据表
),(i i y x i i i y x p r -=)(i i i y x p r -=)(i
m
i r ≤≤0max T
m r r r r ),,(10 =∑
=m
i i
r 0
∑=m
i i
r
2∑=m
i i
r
02
i r
六、程序流程图
七、实验结果 >> zuixiaoerchenfa ans =
27-May-2015 ans =
7.3611e+05 ans =
1.0e+03 *
2.0150 0.0050 0.0270 0.0140 0.0010 0.0213 >>
x
y
x
y
八、实验结果分析 实验程序 quxiannihe.m clear all
date,now,clock
x0=[0.0 0.5 0.6 0.7 0.8 0.9 1.0];
y0=[1 1.75 1.96 2.19 2.44 2.71 3.00];
w=ones(size(x0));
x=0:0.01:1;
%进行五次曲线拟合
N=5;
for i=1:N
a1=LSF(x0,y0,w,i) ;
y=polyval(a1,x);
figure(i)
plot(x0,y0,'ok',x,y,'r')
title('最小二乘法');
最小二乘法
x
y
x
y
legend('y0','y');
xlabel('x');
ylabel('y');
end
实验二：4阶经典龙格库塔法解常微分方程
一、实验名称：4阶经典龙格库塔法解常微分方程
实验时间： 2015-5-14
实验地点：主楼机房
实验器材：计算机matlab
二、实验目的：学习掌握4阶经典R-K方法，体会参数和步长对问题的影响。
三、实验要求：
（1）用4阶经典R-K法编写计算程序，要求用法与ode45一致。
并将计算结果画图比较，并分析步长变化对解的影响。
（2）当激励力幅值F分别按0.3，0.33，0.4， 0.43，0.54，0.58，0.75， 0.84，11.21，13.34进行计算。
每一个数据画出三幅图，分别为时间位移曲线，时间速度曲线和相图。
考察激励力幅值F变化引起的系统响应的变化。
（3）请采用MATLAB中的内部库函数ode45求解此常微分方程初值问题的解，并与（1）中的结果进行比较。
四、算法描述（实验原理与基础理论）
系统方程和表述如下：
则系统的输出按如下求解：
其中：
这样，下一个值(y n +1)由现在的值(y n )加上时间间隔(h )和一个估算的斜率的乘积决定。
该斜率是以下斜率的加权平均：
k 1是时间段开始时的斜率；
k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y 在点tn + h/2的值； k3也是中点的斜率，但是这次采用斜率k2决定y 值； k4是时间段终点的斜率，其y 值用k3决定。
五、实验内容：求解常微分方程初值问题，考虑著名的Duffing 方程。
G ．Duffing 在1918 年引入了一个带有立方项的非线性振子来描述出现在许多力学问题中的质量、弹簧、阻尼系统。
从那时起，Duffing 方程在非线性动力学系统的研究中占有重要的地位。
Duffing 方程的标准形式是
22d d ()()d d x x
c f x g t t t
++= 其中：()f x 是一个含有三次项的非线性函数，()g t 是一个周期函数。
把3()f x x x =-，()cos()g t F t ω=代入上式，可得
223d d cos()d d x x x c F t t t
x ω++=- （1） 式中：0.301, 1.201c ω==, 步长0.01h =; 初值向量为:x0=(0, 0.1)。
要求考察激励力幅值F 变化引起的系统响应的变化。
积分时间区间为：[0, 50]。
六、程序流程图
开始
输入a,b,n,x a,y y0
七、实验结果
i x(i) y(i)
1 0.0000 1.0000
2 0.1000 0.9901
3 0.2000 0.9615
4 0.3000 0.9174
5 0.4000 0.8621
6 0.5000 0.8000
7 0.6000 0.7353
8 0.7000 0.6711
9 0.8000 0.6098
10 0.9000 0.5525
11 1.0000 0.5000
12 1.1000 0.4525
13 1.2000 0.4098
>>
八、实验结果分析
实验程序
function varargout=saxplaxliu(varargin)
clc,clear
x0=0;xn=1.2;y0=1;h=0.1;
[y,x]=lgkt4j(x0,xn,y0,h);
n=length(x);
fprintf('i x(i) y(i)\n');
for i=1:n
fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i) ); end
function z=f(x,y)
z=-2*x*y^2;
function [y,x]=lgkt4j(x0,xn,y0,h)
x=x0:h:xn;
n=length(x);
y1=x;
y1(1)=y0;
for i=1:n-1
K1=f(x(i),y1(i));
K2=f(x(i)+h/2,y1(i)+h/2*K1);
K3=f(x(i)+h/2,y1(i)+h/2*K2);
K4=f(x(i)+h,y1(i)+h*K3);
y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);
end
y=y1;。