数值计算方法第3、4次实验--龙贝格--龙格库塔
龙格库塔解答
实验题目3 四阶Runge-Kutta 方法摘要00(,)()dyf x y dxy x y⎧=⎪⎨⎪=⎩（6.1） 常微分方程初值问题的数值解法是求方程（6.1）的解在点列1(0,1,)n n n x x h n -=+= 上的近似值n y ，这里n h 是1n x -到n x 的步长，一般略去下标记为h 。
经典的R K -方法是一个四阶的方法，它的计算公式是：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （6.2）R K-方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f 。
在用龙格库塔方法时,要注意N 的选择要合适，N 太大，会使计算量加大，N 太小，h 较大，可能会使误差增大。
因此选择合适的N 很重要。
我们要在考虑精度的基础上，选择合适的N 。
在此，用c 语言实现了龙格库塔方法。
前言在此程序中把N 定义成宏常量，运行该程序，只需输入x0，y0，修改函数f ，再运行程序即可。
龙格-库塔法流程图程序设计流程#include<stdio.h>#define N 5double f(double x,double y);void main(){double a=0,b,aa,h,x,y,K1,K2,K3,K4;long n;printf("please input a,b,aa\n");scanf("%lf,%lf,%lf",&a,&b,&aa);h=(double)(b-a)/N;for(n=1;n<=N;n++){K1=h*f(a,aa);K2=h*f(a+h/2,aa+K1/2);K3=h*f(a+h/2,aa+K2/2);K4=h*f(a+h,aa+K3);x=a+h;y=aa+(double)1/6*(K1+2*K2+2*K3+K4);printf("x1=%lf,y1=%lf\n",x,y);a=x;aa=y;}}double f(double x,double y){double f;f=x+y;return f;}龙格库塔思考题解答:问题1解答:数值解和解析解不一样,因为数值解是对给定的x值求出y值,而没有给出y与x的函数关系.相反解析解就给出了y与x的函数关系。
龙格库塔法
一、高阶泰勒法
假设初值问题
龙格—库塔法 龙格 库塔法
dy = f (t , y ) dt y (a) = α 的解y (t)及f (t , y )足够光滑.
将y (ti +1 )在ti处作n阶泰勒展开, 得
a≤t ≤b
(1)
y′′(ti ) 2 y ( n ) (ti ) n y ( n +1) (ξ i ) n +1 y (ti +1 ) = y (ti ) + y′(ti )h + h +L+ h + h n! 2! (n + 1)! 其中, ti < ξ i < ti +1.
2
i
i
1
3
i
i
2
4
i
i
3
i +1
i
6123 Nhomakorabea4
作业 教材P198 习题3
(2)
(3)
首先将y (ti +1 )在ti处展成幂级数 h2 y (ti +1 ) = y (ti ) + hy′(ti ) + y′′(ti ) + O(h 3 ) 2 将 y′(t ) = f (t , y (t )) y′′(t ) = f t′(t , y (t )) + f y (t , y (t )) f (t , y (t )) 代入上式, 得 h2 y (ti +1 ) = y (ti ) + hf + ( f t + ff y ) + O(h 3 ) (3) 2 其中f , f t , f y′分别表示相应函数在点(ti , y (ti ))处的函数值.
龙格库塔
数值分析中，龙格－库塔法（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阶方法控制误差。
龙格库塔法介绍
yn
hf
(xn, yn ))]，
(x, y,h) 1[ f (x, y) f (x h, y hf (x, y))],
2
|
( x,
y1,
h)
(x,
y2 ,
h)
|
[L
2
L 2
(1
hL)]
|
y1
y2
|,
L
L(1
h0L)，h 2
h0.
类似地,不难验证其他龙格 库塔方法的收敛性.
这里c1,c2,c3,2,3, 21, 31, 32均为待定参数.
Tn1 y(xn1) yn1 O(h4 )
(3.11)
c1 c2 c3 1
2
21
3 31 32
c22
c33
1 2
cc232223c2332
将步长折半，从xn用两步求xn1处的近似值，则有
y(xn1)
h
yn21
2c
h 2
5
.
从而
h
y ( xn 1) y ( xn 1)
yn21 ynh1
1， 16
得到事后估计式：
y ( xn 1)
h
yn21
1 15
(
h
yn21
ynh1).
通过检查步长折半前后计算结果的偏差，
y(x) (x, y(x),0) 0 p 1 单步法(4.1)收敛. 定义4 若单步法(4.1)增量函数(x, y,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经过上述修改和改写后，文章已经没有格式错误和明显的问题段落了。
龙格库塔实验报告
一、实验背景常微分方程（ODE）在自然科学、工程技术等领域中具有广泛的应用。
然而，许多微分方程无法得到精确解析解，因此需要借助数值方法进行求解。
龙格-库塔（Runge-Kutta）方法是一种常用的数值求解常微分方程的方法，具有精度高、稳定性好等优点。
本实验旨在通过编写程序，实现四阶龙格-库塔方法，并验证其在求解常微分方程中的有效性和准确性。
二、实验目的1. 理解四阶龙格-库塔方法的基本原理和计算步骤。
2. 编写程序实现四阶龙格-库塔方法。
3. 选取典型常微分方程，验证四阶龙格-库塔方法的求解精度和稳定性。
三、实验原理四阶龙格-库塔方法是一种基于泰勒级数展开的数值方法，其基本思想是将微分方程的解在某个区间内进行近似，并通过迭代计算得到近似解。
具体步骤如下：1. 初始化：给定初始条件y0，步长h，求解区间[a, b]。
2. 迭代计算：对于k=1, 2, ..., n（n为迭代次数），- 计算k1 = f(xk-1, yk-1)（f为微分方程的右端函数）；- 计算k2 = f(xk-1 + h/2, yk-1 + h/2 k1)；- 计算k3 = f(xk-1 + h/2, yk-1 + h/2 k2)；- 计算k4 = f(xk-1 + h, yk-1 + h k3)；- 更新y值：yk = yk-1 + (h/6) (k1 + 2k2 + 2k3 + k4)；- 更新x值：xk = xk-1 + h；3. 输出结果：输出最终的近似解y(n)。
四、实验步骤1. 编写程序实现四阶龙格-库塔方法。
2. 选取典型常微分方程，如：- y' = -y，初始条件y(0) = 1，求解区间[0, 2π]；- y' = y^2，初始条件y(0) = 1，求解区间[0, 1]。
3. 对每个常微分方程，设置不同的步长h和迭代次数n，分别计算近似解y(n)。
4. 将计算得到的近似解与解析解进行比较，分析四阶龙格-库塔方法的精度和稳定性。
河北工业大学数值分析实验三实验四实验报告
数值分析实验报告指导老师：宛艳萍姓名：班级：学号：实验三 复化辛卜生法，龙贝格法1.实验名称：复化辛卜生法，龙贝格法2.实验目的1）通过实际计算体会各种方法的精确度。
2）会编写用复化辛卜生、龙贝格算法求定积分的程序。
3.算法描述1）用复化辛卜生法计算积分 dxx I ⎰+=12)1/(1算法：复化辛卜生公式为S n =h/6∑∑+-=+++)]()2/(4)([11k k kn k x f h x f xf ,计算过程为：1．令,/)(n a b h -= ),2/(1h a f s +=;02=s2．对1,,2,1-=n k计算),2/(11h kh a f s s +++=)(22kh a f s s ++=3．))(24)((6/21b f s s a f h s +++= 。
2）龙贝格算法计算dxxI ⎰+=102)1/(156e ε=-算法)((12/12∑-=++=n k k n n n x f h T T ;/)(n a b h n -= n k h k x )2/1(2/1+=+)(3/122n n n n T T T S -+= )_(15/122n n n n S S S C +=)(63/122n n n n C C C R -+=用事后估计法控制精度2|5e -6n n R R -< 。
4.源程序：1）/* 用复化辛卜生公式求积分 */ #include "stdio.h" float fx(float x){double f;f=1.0/(1.0+x*x); return f; } double fs(int n){double a=0.0,b=1.0,h,s,s1,s2=0; int i;h=(b-a)/n; s1=fx(a+h/2); for(i=1;i<n;i++){s1=s1+fx(a+i*h+h/2); s2=s2+fx(a+i*h);}s=(h/6.0)*(fx(a)+fx(b)+4*s1+2*s2);return s;}void main(){printf("实验三复化辛卜生法计算机112 耿向飞学号：112434\n");printf("s(2)=%lf\ns(4)=%lf\ns(8)= %lf",fs(2),fs(4),fs(8));}2）/* 龙贝格法 */#include "stdio.h"#include "math.h"#define E 2.71828182//被积函数f(x)double fx(double x){double f;f=1/(1+x*x);return f;}//梯形公式求tndouble tx(int n){double s3=0.0,h,t,b=1.0,a=0.0;int i;h=(b-a)/n;for(i=1;i<n;i++)s3=s3+fx(i*h);t=(h/2)*(fx(a)+fx(b)+2*s3);return t;} double s(int n){double s;s=tx(2*n)+(1.0/3.0)*(tx(2*n)-tx(n ));return s;}double c(int n){double c;c=s(2*n)+(1.0/15.0)*(s(2*n)-s(n)) ;return c;}double r(int n){double r;r=c(2*n)+(1.0/63.0)*(c(2*n)-c(n)) ;return r;}void main(){double rr,pp;int n=1;rr=r(n);pp=r(2*n)-r(n);printf("实验三龙贝格法计算机112 耿向飞学号：112434\n");printf("结果为：%.15lf 误差小于等于： %.15lf",rr,pp);}5.运行结果1）复化辛卜生公式2）龙贝格算法6.对算法的理解与分析：复化辛卜生公式和龙贝格算法适用于求数值积分，而且都能提高计算积分的精度龙贝格算法其实是在复化辛卜生公式递推的基础之上生成的一种精度高，而且收敛速度也较快的一种算法。
实验龙格—库塔方法
实验 龙格—库塔方法实验目的： 深入理解龙格—库塔方法的原理，同时学会用迭代方法看待某些特定问题。
比较本方法与其它方法（尤其是欧拉方法）解题的不同处，分析两者的精度。
实验要求: 编写程序，用龙格—库塔方法求解第三章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 阶精度，只需
龙格-库塔方法(Runge-Kutta)
龙格-库塔方法〔Runge-Kutta〕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值(即由前面(i-1)个已算出的表示，故公式是显式的.例)，而ki如当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 。
龙格-库塔方法
h 按照数学学 习一般方法， yn 1 yn （k1 2k 2 2k3 +k 4） 6 继续上面的 做法，经过 k1 f ( xn , yn ) 较复杂的数 h 学推理和计 k （16） 2 f ( xn 1 , yn k1 ) 算，可以导 2 2 出四阶龙格 h －库塔格式， k f ( x 1 , yn k 2 ) 这里有一个 3 n 2 2 较有用的四 k f ( x , y hk ) 阶格式。 n 1 n 3 4
？
3、欧拉格式
yn 1 yn hf ( xn , yn ),
取点 xn 的斜率 k
*
(3)
n 0,1, 2,
f ( xn , yn )作为区间 [ xn , xn1 ] 上的平均斜率，精度低。
4、改进的欧拉格式
于是：
h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn hf ( xn , yn ))] 2 h yn 1 yn [k1 k2 ] 2
龙格－库塔方法
考虑一维经典初值问题
dy f ( x , y ) , y( x0 ) y0 , x [a, b] dx
(1)
龙格－库塔方法
一、设计思想
加权平均斜率。
1、微分中值定理
如果 f (x) 在 [a, b]上连续，在 ( a, b) 内可导，则在( a, b) 内至 少存在一点 ，有 f (b) f (a) f ' ( ) ba
xnq xn qh, 0 p q 1
令
yn1 yn h[(1 )k1 k2 k3 ] (14)
k1 , k 2
计算方法 实验报告 拉格朗日 龙贝格 龙格库塔
主界面：
/*lagrange.c*/
float real_value(float x) /*由被插值函数计算真实值*/
c=getchar();
if(c=='N'||c=='n') break;
}
}
/*romberg.c*/
double function(double x) /*被积函数*/
{
return 4.0/(1+(x)*(x));
}
double t(double a,double b,int m) /*计算T1*/
实验二（龙贝格公式）
§公式
§算法描述
§流程图
§运行结果
§结果分析：Romberg积分法是在积分步长逐步折半的过程中，用低精度求积公式的组合得到更高精度求积公式的一种方法，它算法简单，且收敛加速效果极其显著。
实验三（四阶龙格库塔）
§公式
k1=h*f(xn,yn);
k2=h*f(xn+h/2,yn+k1/2);
T1=t(a,b,0);
T2=T1/2.0+t(a,b,1);
S1=(4*T2-T1)/3.0;
T1=T2;
T2=T1/2.0+t(a,b,2);
S2=(4*T2-T1)/3.0;
C1=(16*S2-S1)/15.0;
T1=T2;
T2=T1/2.0+t(a,b,3);
S1=S2;
S2=(4*T2-T1)/3.0;
四阶龙格-库塔法求解常微分方程的初值问题
实验3 龙贝格积分1. 算法原理龙贝格积分法是将积分区间[a, b]逐次分半，利用复化梯形求积公式逐次递推计算，再利用复化辛普森求积公式、复化科茨求积分公式、龙贝格积分公式逐次计算，以得到较高的精度的积分近似值。
计算公式如下：1[()()]2b a T f a f b -=+， 12112211((21))222kk k k k i b a b a T T f a i +++=--=++-∑，0,1,2k =⋅⋅⋅， 111111222222222322221()411()411()41k k k k k k k k k k k k S T T T C S S S R C C C ++++++⎧=--⎪-⎪⎪=--⎨-⎪⎪=--⎪-⎩, 0,1,2k =⋅⋅⋅，计算过程按右图所示的顺序进行，若122||k k R R ε+-<，则取12[]k I f R +=；否则，继续计算，直到满足精度要求为止。
2. 程序框图3.程序使用说明本程序使用MATLAB来求解龙贝格积分。
源程序文件“romber.m”为龙贝格积分的源程序，tscr为T、S、C、R数表，q 为积分值，err为积分误差，f为被积函数，a为积分上限，b为积分下限，ep为精度要求。
输入被积函数f，积分上限a，积分下限b，精度要求ep后，在命令窗口输入[tscr,q,err] = romber(f,a,b,ep)，回车后即可算出T、S、C、R数表、积分值和积分误差。
源程序文件“romber2.m”是计算实习6.2算例的程序，直接运行后，在命令窗口输入[tscr,q,err] = romber(f,a,b,ep)，回车后即可得到算例结果。
4.算例计算结果此算例为课本211页计算实习6.2.首先，设定MATLAB的数值格式为“long”。
假设精度要求ep = 10-13，运行后所得结果如下：tscr =0.750000000000000 0 0 00.708333333333333 0.694444444444444 0 00.697023809523809 0.693253968253968 0.693174603174603 00.694121850371850 0.693154530654531 0.693147901481235 0.6931474776448320.693391202207527 0.693147652819419 0.693147194297078 0.6931471830719330.693208208269249 0.693147210289823 0.693147180787850 0.6931471805734180.693162438883403 0.693147182421455 0.693147180563564 0.6931471805600040.693150995228108 0.693147180676343 0.693147180560002 0.693147180559945 q =0.693147180559945err =5.828670879282072e-14所以积分值为0.693147180559945，误差为5.828670879282072 10-14。
三阶_龙格-库塔公式 详细推导过程
作者最近在学习《数值分析》这门课程时，发现多数教材只以二阶二段龙格库塔（R-K ）公式为例进行推导，对于三四阶R-K 公式只是简单一提，便给出结论。
作者查了一些资料，都未找到三阶R-K 公式的详细推导过程。
于是饶有兴致地尝试推导，功夫不负有心人，终于有所收获，特与大家分享。
作者以30岁的年龄、18岁的心态，终得此成果，很欣慰。
R-K 公式应用广泛，是一种高精度的单步法，读者在阅读本文之前，应该已经了解R-K 方法，并已熟悉二阶二段龙格库塔（R-K ）公式的推导过程。
主要参考资料：同济大学出版社：《高等数学》、《数值分析基础》。
三阶R-K 公式如下：k y iri in∑=++=11n yω 其中：f kn 1),(h hf y x nn ==),(1212k k b y a x nn hf ++= ｈ２ ),(23213133kk b k b y a x nn hf +++= ｈ将k 2按二元函数在处按Taylor 公式展开， （),(f n y x nn f =）)(221),(3222212122222212n121h ff hb ff b ha fh a ff hb fha fk b y a x O n n f yynxyxxy x nn ++++++=++）（ ｈ２然后k 3作类似的处理(注意：将其中的k 2用上式代替，并注意略去导致最终展开式中的O(h 4)项)：)()2)()222122),(322323122232231322331232232132y32313n232313)(h ff hb b ff h b b f b ha ff b ha fh a ff f h b b f f h a f h b f f hb fha fk b f h b y a x O n n n x n n n f yynyynxyxxy yyxnn ++++++++++++=++++（ ｈ将上述公式带入总公式即可，读者可自己尝试，此处不再赘述。
四阶龙格库塔实验报告
三、四阶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 。
龙格库塔实验报告
龙格库塔实验报告龙格库塔实验报告引言：龙格库塔法（Runge-Kutta method）是一种常用的数值求解常微分方程（ODE）的方法。
它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）于19世纪末独立发展而来。
龙格库塔法通过将微分方程转化为一系列逼近值，从而实现对微分方程的数值求解。
本实验旨在通过对龙格库塔法的研究和实践，深入了解其原理和应用。
一、龙格库塔法的原理龙格库塔法的基本思想是将微分方程的解分解为一系列逼近值，然后通过逐步迭代的方式计算这些逼近值。
具体而言，龙格库塔法将求解微分方程的过程分为多个步骤，每个步骤都计算一个逼近值，并利用这些逼近值逐步逼近真实解。
二、龙格库塔法的步骤1. 确定微分方程和初始条件：首先，需要明确待求解的微分方程以及初始条件。
微分方程可以是一阶或高阶的，初始条件是方程在某一点上的已知值。
2. 确定步长：步长（或称为时间间隔）决定了逼近值的精度。
较小的步长能够提高逼近值的准确性，但也会增加计算量。
因此，需要根据具体问题的需求选择合适的步长。
3. 迭代计算：根据龙格库塔法的公式，进行逐步的迭代计算。
首先，根据初始条件计算出第一个逼近值。
然后，利用该逼近值和微分方程的导数计算出下一个逼近值。
重复这一过程，直到达到所需的迭代次数或满足精度要求。
4. 计算误差：为了评估逼近值的准确性，需要计算误差。
误差可以通过与解析解的比较得出，或者通过比较不同步长下的逼近值得出。
较小的误差表明逼近值较为准确。
三、龙格库塔法的应用龙格库塔法广泛应用于科学和工程领域，特别是在求解微分方程模型的数值计算中。
以下是一些常见的应用场景：1. 物理学：龙格库塔法可以用于模拟物体在重力场中的运动，如自由落体、抛体运动等。
通过求解微分方程，可以得到物体的位置、速度等随时间的变化规律。
2. 经济学：龙格库塔法可以用于经济学模型的求解，如经济增长模型、投资模型等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1.完成复合梯形、复合辛普森求积公式，用变步长、事后误差估计的方法完成P145例题6.3.1对应表6-3的计算，课后作业题7.%复合梯形求积公式的MATLAB实现%姓名：王定%学号：1306034248%中北大学仪器与电子学院function T_n=ComTrapezoidal(a,b,n)%a为积分上限%b为积分下限%n为划分区间的个数h=(b-a)/n;for(k=0:n)x(k+1)=a+k*h;if(x(k+1)==0)x(k+1)=10^(-10);endendI_1=h/2*(f(x(1))+f(x(n+1)));for(i=2:n)F(i)=h*f(x(i));endI_2=sum(F);I_n=I_1+I_2%****************************% %复合辛普森求积公式的MATLAB实现function I=ComSimpson(a,b,n)n=2;h=(b-a)/2;I1=0;I2=(f(a)+f(b))/h;eps=1.0e-4;while(abs(I2-I1)>eps)n=n+1;h=(b-a)/n;I1=I2;I2=0;for(i=0:n-1)x=a+h*i;x1=x+h;I2=I2+h/6*(f(x)+4*f((x+x1)/2)+f(x1));endendI=I2 %变步长梯形求积法的MATLAB实现function AdaptiveTrapezoidal(a,b,eps)m=1t0=0;t2=(f(a)+f(b))/4+f((a+b)/2)while(abs(t2-t0)>eps)m=m+1p=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=t1t0=p;endI=t2%**************************%>>AdaptiveTrapezoidal(1,2,0.000001)m = 1t2 = 3.03948481584447m = 2t2 = 2.02304986763725m = 3t2 = 2.02080858246806m = 4t2 = 2.02024624995310m = 5t2 = 2.02010553934348m = 6t2 = 2.02007035369492m = 7t2 = 2.02006155678257m = 8t2 = 2.02005935752321m = 9t2 = 2.02005880770641I = 2.02005880770641%课后作业题7.>>AdaptiveTrapezoidal(1,3,0.0001)m = 1t2 =7.99930630234471m = 2t2 =10.84204346755743m = 3t2 =10.92309388961378m = 4t2 =10.94339842118680m = 5t2 =10.94847716722413m = 6t2 =10.94974701694155m = 7t2 =10.95006448956964m = 8t2 =10.95014385836406I =10.950143858364062.完成龙贝格积分，计算P150例题6.4.2对应的表6-5，按例题列表格显示数值，完成课后作业题8.%龙贝格求积法的MATLAB实现%姓名：王定%学号：1306034248%中北大学仪器与电子学院function [I,step]=Romberg(f,a,b,eps)f=input('please input f=');a=input('please input a=');b=input('please input b=');eps=input('please input eps=');%I为所求积分值%step为积分划分的子区间次数%f为被积函数%a为积分上限%b为积分下限%eps为积分精度if(nargin==3)eps=1.0e-4;end;M=1;tol=10;k=0;T=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsy m(sym(f)),b));while(tol>eps)k=k+1;h=h/2;Q=0;for(i=1:M)x=a+h*(2*i-1);Q=Q+subs(sym(f),findsym(sym(f)),x);endT(k+1,1)=T(k,1)/2+h*Q;M=2*M;for(j=1:k)T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1);endtol=abs(T(k+1,j+1)-T(k,j));endI=T(k+1,k+1)step=k %计算P150例题6.4.2对应的表6-5>> [I,step]=Romberg('4/(1+x^2)',0,1); format longI,stepI =3.14159266527772step =4%********************************%%按例题列表格显示数值，完成课后作业题8.>>please input f='(2/sqrt(pi))*exp(-x)'please input a=0please input b=1please input eps=1.0e-5I =0.71327166981418step =31、完成预报校正公式的程序，完成P183例题7.3.2对应表7—6第3列结果的计算，课后作业题3。
%预报校正公式的MATLAB实现%姓名：王定%学号：1306034248%中北大学仪器与电子学院functiony=predict_collect(f,h,a,b,y0,varvec)%f为一阶常微分方程一般表达式的右端项%h为积分步长%a为自变量的取值下限%b为自变量的取值上限%y0为函数初值%varvec为常微分方程变量组clearclcsyms x y;z=y-2*x/y;f=input('please input f=');h=input('please input h=');a=input('please input a=');b=input('please input b=');y0=input('please input y0=');varvec=input('please input varvec=');%you can input [x y]varvec=[x y];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)+hy(i-1)+h*K1]);y(i)=y(i-1)+h*(K1/2+K2/2);end %完成P183例题7.3.2对应表7—6第3列结果的计算please input f=zplease input h=0.1please input a=0please input b=1.4please input y0=1please input varvec=[x y]ans =1.000000000000001.095909090909091.184096569243001.266201360875781.343360151484001.416401928536911.485955602415671.552514091326151.616474782752061.678166363675191.737867401035411.795819744910661.852238599050291.90732041783776%课后作业题3please input f=x+yplease input h=0.1please input a=0please input b=1please input y0=1please input varvec=[x y]ans =1.000000000000001.110000000000001.242050000000001.398465250000001.581804101250001.794893531881252.040857352728782.323147374765302.645577849115663.012363523272813.428161693216452、完成四阶标准龙格—库塔公式的程序，完成P183例题7.3.2对应表7—6第4、5列结果的计算，课后作业题第4、6题。
%四阶龙格-库塔法的MATLAB实现%姓名：王定%学号：1306034248%中北大学仪器与电子学院function y=Runge_Kutta4(f,h,a,b,y0,varvec) %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/2y(i-1)+K1*h/2]);K3=Funval(f,varvec,[x(i-1)+h/2y(i-1)+K2*h/2]);K4=Funval(f,varvec,[x(i-1)+hy(i-1)+h*K3]);y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6; endformat short; %P183例题7.3.2对应表7—6第4列结果的计算>> syms x y;z=y-2*x/y;result=Runge_Kutta4(z,0.2,0,0.6,1,[x y]);>> format longresult =1.000000000000001.183229287445311.34166692985261%*********************%%P183例题7.3.2对应表7—6第5列结果的计算>> syms x y;z=y-2*x/y;result=Runge_Kutta4(z,0.05,0,0.6,1,[x y]);>> format longresult =1.000000000000001.048808861191141.095445139884111.140175461296841.183216003999081.224744930067011.264911134342841.303840563372121.341640881476441.378404983543001.414213684897401.44913781228807%课后作业题第4题>> syms x y;z=x+y;result=Runge_Kutta4(z,0.2,0,1,1,[x y]);>> format long>> resultresult =1.000000000000001.242800000000001.583635920000002.044212912688002.651041651557123.43650227321187%*********************%%课后作业题第6题>> syms x y;z=x^2-y^2;result=Runge_Kutta4(z,0.1,-1,0,0,[x y]);format longresultresult =0.090047357462400.160726883901280.213482650382680.250368369544590.273775247604510.286221918318140.290213341573770.288160170936890.282343663961140.274910519234623、用ode23和ode45计算1、2中的题目%用ode23完成P183例题7.3.2对应表7—6第3、4、5列结果的计算>>[t,y]=ode23(@func,[0,1.4] ,1)t =0.080000000000000.220000000000000.360000000000000.500000000000000.640000000000000.780000000000000.920000000000001.060000000000001.200000000000001.340000000000001.40000000000000y =1.000000000000001.077036361346181.200031136790371.311546822967341.414303716907341.510093464126161.600170812442271.685455644781501.766646865702611.844291109230861.918826503754141.94990710209889 %用ode23完成课后作业题3>>[t,y]=ode23(@func,[0,1],1)t =0.080000000000000.180000000000000.280000000000000.380000000000000.480000000000000.580000000000000.680000000000000.780000000000000.880000000000000.980000000000001.00000000000000y =1.000000000000001.086570666666671.214421681777781.366235028644741.544530745823881.752093895926361.992002437314612.267658026972202.582820062808782.941643306080833.348719460437003.43636669849711%用ode23完成P183例题7.3.2对应表7—6第3、4、5列结果的计算%用ode23完成课后作业题3%用ode45完成P183例题7.3.2对应表7—6第3、4、5列%用ode45完成课后作业题3结果的计算%用ode45完成P183例题7.3.2对应表7—6第3、4、5列结果的计算>> [t,y]=ode45(@func,[0,0.6],1) %用ode45完成课后作业题3 >>[t,y]=ode45(@func,[0,1],1)t =0.035000000000000.070000000000000.105000000000000.140000000000000.175000000000000.210000000000000.245000000000000.280000000000000.315000000000000.350000000000000.385000000000000.420000000000000.455000000000000.490000000000000.525000000000000.560000000000000.595000000000000.630000000000000.665000000000000.700000000000000.735000000000000.770000000000000.805000000000000.840000000000000.875000000000000.910000000000000.945000000000000.980000000000001.015000000000001.050000000000001.085000000000001.120000000000001.155000000000001.190000000000001.225000000000001.260000000000001.295000000000001.330000000000001.365000000000001.40000000000000 y =1.000000000000001.034408313698181.067707873647801.099999811523621.131370863523241.161895104725561.191637539132891.220655493598291.248999620145811.276714583380471.303840484014931.330413445440531.356466022151961.382027529581721.407124732650101.431782103619121.456022008361561.479864888490641.503329647925861.526433763432251.549193374741161.571623396650611.593737762729231.615549464951841.637070597777071.658312433330501.679285589651481.700000034539531.720465105937351.740689565923681.760681725367921.780449429267861.800000064303671.819340599782731.838477685157591.857417625820221.876166383472341.894729608547341.913112719737051.931320875080911.94935896824202t =0.025000000000000.050000000000000.075000000000000.100000000000000.125000000000000.150000000000000.175000000000000.200000000000000.225000000000000.250000000000000.275000000000000.300000000000000.325000000000000.350000000000000.375000000000000.400000000000000.425000000000000.450000000000000.475000000000000.500000000000000.525000000000000.550000000000000.575000000000000.600000000000000.625000000000000.650000000000000.675000000000000.700000000000000.725000000000000.750000000000000.775000000000000.800000000000000.825000000000000.850000000000000.875000000000000.900000000000000.925000000000000.950000000000000.975000000000001.00000000000000y =1.000000000000001.025630247832681.052542202197921.080768305822271.110341836666671.141296914159361.173668496437661.207492438259491.242805517459491.279645441837501.318050846110251.358061355913461.399717617040431.443061303054941.488135111922351.534982836742891.583649398065261.634180852300191.686624387997301.741028404017381.797442545244481.855917711884351.916506055349711.979261064653972.044237605879242.111491932452662.181081680598902.253065964820042.327505421514432.404462220333802.484000059157502.566184269614072.651081865287582.738761554269112.829293733602482.922750605903613.019206232636723.118736547985963.221419352716863.327334443060803.43656366959418。