四阶龙格库塔法解微分方程

合集下载

c++经典四阶龙格库塔法解一阶微分方程组

首先，龙格库塔法（Runge-Kutta method）是一种常用的数值积分算法，它可以用来数值地求解一阶常微分方程（ODE）。

对于一维ODE，经典的四阶龙格库塔法（RK4）是最常用的方法之一。

下面是一个使用C++实现经典四阶龙格库塔法解一阶微分方程组的示例代码：cpp#include <iostream>#include <vector>// 定义微分方程组dy/dx = f(x, y)std::vector<double> f(double x, const std::vector<double>& y) {std::vector<double> dy(y.size());// 此处为具体的微分方程组表达式，请根据实际需求来定义dy[0] = x * y[0]; // 示例中假设有一个一维微分方程y' = x*yreturn dy;}// 经典四阶龙格库塔法void rk4(double x0, double y0, double h, int numSteps) {double x = x0;double y = y0;for (int i = 0; i < numSteps; ++i) {std::vector<double> k1 = f(x, std::vector<double>{y});std::vector<double> k2 = f(x + h / 2, std::vector<double>{y + h * k1[0] / 2});std::vector<double> k3 = f(x + h / 2, std::vector<double>{y + h * k2[0] / 2});std::vector<double> k4 = f(x + h, std::vector<double>{y + h * k3[0]});y += h * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6;x += h;std::cout << "x: " << x << ", y: " << y << std::endl;}}int main() {double x0 = 0.0; // 初值xdouble y0 = 1.0; // 初值ydouble h = 0.1; // 步长int numSteps = 10; // 迭代次数rk4(x0, y0, h, numSteps);return 0;}请注意，这只是一个简单的示例代码，用于演示如何使用经典四阶龙格库塔法求解一维微分方程组。

四阶龙格库塔实验报告

三、四阶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 。

四阶龙格库塔法例题

四阶龙格库塔法例题四阶龙格库塔法是一个常见的求解常微分方程（ODE）的数值方法。

它具有高精度和稳定性的特点，常被应用于物理、化学、生物等领域中的ODE求解问题。

下面我们来看一个四阶龙格库塔法的例题。

首先，我们需要了解ODE的定义和求解过程。

ODE是描述自变量和它的某个函数之间关系的一个方程，其中自变量通常是时间。

求解ODE通常通过数值方法，将连续的函数曲线离散化成若干个点，然后计算它们之间的差值和斜率，最终得到一条近似的曲线。

四阶龙格库塔法是一种基于差积公式的ODE求解数值方法，通过多次迭代得到数值解的近似值。

接下来，我们看一个例子：已知ODE y' = y + x，初始值 y0=0，则求y(0.1)的近似数值解。

我们可以将其转化为差分方程：y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6,其中h=0.1，k1=f(i,y(i)),k2=f(i+h/2,y(i)+k1*h/2), k3=f(i+h/2,y(i)+k2*h/2),k4=f(i+h,y(i)+k3*h)，其中f(i,y(i))为ODE的右侧，即y(i)+x，代入计算得到：k1=0+0=0，k2=0.055，k3=0.057，k4=0.116，因此y(0.1)=0+0.05*(0+2*0.055+2*0.057+0.116)/6=0.0058。

以上就是四阶龙格库塔法例题的具体过程。

需要注意的是，在实际应用中，我们需要根据ODE的具体情况确定差分方程的形式以及步长h的大小，以保证求解的精度和稳定性。

此外，还需要注意初始值和边界条件的设定，以避免数值误差的累积和算法发散等问题的出现。

总之，四阶龙格库塔法是数值计算中基本的数值方法之一，具有广泛的应用和重要的意义。

掌握其基本理论和具体运用，可以帮助我们更好地解决现实生活和科学研究中的相关问题。

四阶龙格-库塔(R-K)方法求常微分方程

中南大学MATLAB程序设计实践材料科学与工程学院2013年3月26日一、编程实现“四阶龙格－库塔（R-K）方法求常微分方程”，并举一例应用之。

【实例】采用龙格-库塔法求微分方程：⎩⎨⎧==+-=0, 0)(1'00x x y y y 1、算法说明：在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。

其算法公式为：)22(63211k k k hy y n n +++=+其中：⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++=++==) ,()21,21()21 ,21(),(3423121hk y h x f k hk y h x f k hk y h x f k y x f k n n n n n n n n 2、流程图：2.1、四阶龙格－库塔（R-K ）方法流程图：2.2、实例求解流程图：3、源程序代码3.1、四阶龙格－库塔（R-K）方法源程序：function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x))%此程序可解高阶的微分方程。

只要将其形式写为上述微分方程的向量形式%函数 f(x,y): fun%自变量的初值和终值：x0, xt%y0表示函数在x0处的值，输入初值为列向量形式%自变量在[x0,xt]上取的点数：PointNum%varargin为可输入项，可传适当参数给函数f(x,y)%x：所取的点的x值%y：对应点上的函数值if nargin<4 | PointNum<=0PointNum=100;endif nargin<3y0=0;endy(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长x=x0+[0:(PointNum-1)]'*h; %得x向量值for k=1:(PointNum) %迭代计算f1=h*feval(fun,x(k),y(k,:),varargin{:});f1=f1(:)'; %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:});f2=f2(:)'; %得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:});f3=f3(:)'; %得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});f4=f4(:)'; %得公式k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1) end3.2、实例求解源程序：%运行四阶R-K法clear, clc %清除内存中的变量x0=0;xt=2;Num=100;h=(xt-x0)/(Num-1);x=x0+[0:Num]*h;a=1;yt=1-exp(-a*x); %真值解fun=inline('-y+1','x','y'); %用inline构造函数f(x,y)y0=0; %设定函数初值PointNum=5; %设定取点数[x1,y1]=ode23(fun,[0,2],0);[xr,yr]=MyRunge_Kutta(fun,x0,xt,y0,PointNum);MyRunge_Kutta_x=xr'MyRunge_Kutta_y=yr'plot(x,yt,'k',x1,y1,'b--',xr,yr,'r-')legend('真值','ode23','Rung-Kutta法解',2)hold onplot(x1,y1,'bo',xr,yr,'r*')4、程序运行结果：MyRunge_Kutta_x =0 0.5000 1.0000 1.5000 2.0000MyRunge_Kutta_y =0 0.3932 0.6318 0.7766 0.8645二、变成解决以下科学计算问题：（一）[例7-2-4] 材料力学复杂应力状态的分析——Moore 圆。

四阶龙格库塔法解微分方程

四阶龙格库塔法解微分方程一、四阶龙格库塔法解一阶微分方程①一阶微分方程：cos y t ，初始值y(0)=0，求解区间[0 10]。

MATLAB 程序：%%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程%%%%%%%%%%% y'=cost%%%%%%%%%%% y(0)=0, 0≤t ≤10，h=0.01%%%%%%%%%%% y=sinth=0.01;hf=10;t=0:h:hf;y=zeros(1,length(t));y(1)=0;F=@(t,y)(cos(t));for i=1:(length(t)-1)k1=F(t(i),y(i));k2=F(t(i)+h/2,y(i)+k1*h/2);k3=F(t(i)+h/2,y(i)+k2*h/2);k4=F(t(i)+h,y(i)+k3*h);y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h;endsubplot(211)plot(t,y,'-')xlabel('t');ylabel('y');title('Approximation');span=[0,10];p=y(1);[t1,y1]=ode45(F,span,p);subplot(212)plot(t1,y1)xlabel('t');ylabel('y');title('Exact');图1②一阶微分方程：()22*/x t x x t =- ，初始值x(1)=2，求解区间[1 3]。

MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解微分方程%%%%%%%%%%% x'(t)=(t*x-x^2)/t^2%%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128%%%%%%%%%%% 精确解：x(t)=t/(0.5+lnt)h=1/128; %%%%% 步长tf=3;t=1:h:tf;x=zeros(1,length(t));x(1)=2; %%%%% 初始值F_tx=@(t,x)(t.*x-x.^2)./t.^2;for i=1:(length(t)-1)k_1=F_tx(t(i),x(i));k_2=F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1);k_3=F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2));k_4=F_tx((t(i)+h),(x(i)+k_3*h));x(i+1)=x(i)+(1/6)*(k_1+2*k_2+2*k_3+k_4)*h; endsubplot(211)plot(t,x,'-');xlabel('t');ylabel('x');legend('Approximation');%%%%%%%%%%%%%%%%%%%%%%%%%%%% ode45求精确解t0=t(1);x0=x(1);xspan=[t0 tf];[x_ode45,y_ode45]=ode45(F_tx,xspan,x0);subplot(212)plot(x_ode45,y_ode45,'--');xlabel('t');ylabel('x');legend('Exact');图2二、四阶龙格库塔法解二阶微分方程①二阶微分方程：cos y t ，初始值y(0)=0，y'(0)=-1，求解区间[0 10]。

四阶龙格库塔实验报告

4阶经典龙格库塔公式求解微分方程

4阶经典龙格库塔公式求解微分方程4阶龙格库塔法（也称为四阶Runge-Kutta方法）是一种常用于求解常微分方程的数值方法。

它是由德国数学家卡尔·龙格以及他的学生马丁·威尔海姆·库塔于1901年独立提出的。

以下是详细的介绍。

1.问题描述我们考虑一个典型的常微分方程初值问题：dy/dx = f(x, y)，y(x0) = y0其中，x0是给定的初始点，y(x)是我们要求解的函数，f(x,y)是已知的函数。

2.方法原理四阶龙格库塔方法的基本思想是通过使用全区间的函数信息来估计下一步的函数值。

在每个步骤中，我们计算四个不同的估计量，然后将它们组合起来以获取最终的解。

具体来说，我们首先计算在初始点x0上的斜率k1=f(x0,y0)。

然后我们计算在x0+h/2处的斜率k2=f(x0+h/2,y0+h/2*k1)，其中h是步长。

以此类推，我们计算k3和k4分别在x0+h/2和x0+h处的斜率。

然后，我们通过加权组合这些斜率来计算下一个点的函数值y1：y1=y0+(h/6)*(k1+2*k2+2*k3+k4)。

3.算法步骤以下是使用四阶龙格库塔法求解常微分方程的详细步骤：步骤1：给定初始条件 y(x0) = y0，选择步长 h 和积分终点 x_end。

步骤2：计算积分步数n：n = (x_end - x0)/h。

步骤3：设置变量x=x0和y=y0，并将步长设置为h。

步骤4：对每个步数i=1,2,...,n，执行以下计算：4.1：计算斜率k1=f(x,y)。

4.2：计算斜率k2=f(x+h/2,y+h/2*k1)。

4.3：计算斜率k3=f(x+h/2,y+h/2*k2)。

4.4：计算斜率k4=f(x+h,y+h*k3)。

4.5：计算下一个点的函数值y1=y+(h/6)*(k1+2*k2+2*k3+k4)。

4.6：将x更新为x+h，y更新为y1步骤5：重复步骤4，直到达到积分终点 x_end。

常微分方程组的四阶RUNGEKUTTA龙格库塔法MATLAB实现

常微分方程组的四阶RUNGEKUTTA龙格库塔法MATLAB实现欢迎使用 Python 版的实现！数值解常微分方程组的四阶 Runge-Kutta 方法（也被称为 RK4 方法）可以通过迭代的方式逐步逼近精确解。

对于一阶常微分方程组：dy/dt = f(t, y)我们可以通过下面的公式来计算下一个时间步长上的近似解：y(n+1)=y(n)+(1/6)(k1+2k2+2k3+k4)其中k1=h*f(t(n),y(n))k2=h*f(t(n)+h/2,y(n)+k1/2)k3=h*f(t(n)+h/2,y(n)+k2/2)k4=h*f(t(n)+h,y(n)+k3)这里，h代表时间步长，t(n)代表当前时间，y(n)代表当前解。

f(t,y)是给定的方程组。

对于四阶 Runge-Kutta 方法的 MATLAB 实现，可以按照以下步骤进行。

1.首先，定义需要求解的常微分方程组。

function dydt = equations(t, y)dydt = zeros(2, 1);dydt(1) = y(2); % 根据方程组的具体形式修改dydt(2) = -y(1); % 根据方程组的具体形式修改end2.定义RK4方法的求解函数。

function [t, y] = rk4_solver(equations, tspan, y0, h) t = tspan(1):h:tspan(2);y = zeros(length(y0), length(t));y(:,1)=y0;for i = 1:length(t)-1k1 = h * equations(t(i), y(:, i));k2 = h * equations(t(i) + h/2, y(:, i) + k1/2);k3 = h * equations(t(i) + h/2, y(:, i) + k2/2);k4 = h * equations(t(i) + h, y(:, i) + k3);y(:,i+1)=y(:,i)+(k1+2*k2+2*k3+k4)/6;endend3. 调用 rk4_solver 函数求解常微分方程组。

经典四阶龙格库塔法解一阶微分方程组

1.经典四阶龙格库塔法解一阶微分方程组1.1运用四阶龙格库塔法解一阶微分方程组算法分析（1-1），（1-2）（1-3）（1-4）（1-5）（1-6）（1-7）（1-8）（1-9）（1-10）经过循环计算由推得……每个龙格-库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为,一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。

4阶龙格-库塔方法(RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。

1.2经典四阶龙格库塔法解一阶微分方程流程图图1-1 经典四阶龙格库塔法解一阶微分方程流程图1.3经典四阶龙格库塔法解一阶微分方程程序代码：#include <iostream>#include <iomanip>using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h){double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=initial[0];x0=initial[1];y0=initial[2];f1=f(t0,x0,y0); g1=g(t0,x0,y0);f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2); g2=g(t0+h/2, x0+h*f1/2,y0+h*g1/2);f3=f(t0+h/2, x0+h*f2/2,y0+h*g2/2); g3=g(t0+h/2, x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3); x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6; resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<"输入所求微分方程组的初值t0,x0,y0:";cin>>initial[0]>>initial[1]>>initial[2];cout<<"输入所求微分方程组的微分区间[a,b]:";cin>>a>>b;cout<<"输入所求微分方程组所分解子区间的个数step:";cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision(10); H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;dx=x+2*y;return(dx);}double g(double t,double x, double y){double dy;dy=3*x+2*y;return(dy);}1.4经典四阶龙格库塔法解一阶微分方程程序调试结果图示：应用所编写程序计算所给例题：其中初值为求解区间为[0,0.2]。

四阶Runge-Kutta法求解常微分方程

call update_half_gene(p%value,rt_current+rt_step)
call update_fish_state(p%value,rt_current+rt_step)
p=>p%next
end do
!
END SUBROUTINE runge_kutta
解答2：最常用的是四步的Runge-Kutta格式
最佳答案
没试过matlab，算这玩意太慢了，有fortran版的要不，有兴趣的话可以参考一下。
代码:
SUBROUTINE runge_kutta()
!关于Runge-Kutta方法，该方法是用来解形如y'=f(t,y)的常微分方程的
!经典的4阶R-K方法的公式如下：
! Yn+1 = Yn + h/6 * (K1+2K2+2K3+K4)
fk3=f(x(i-1)+h/2,y(i-1)+fk2*h/2);
fk4=f(x(i-1)+h,y(i-1)+fk3*h);
y(i)=y(i-1)+h*(fk1+2*fk2+2*fk3+fk4)/6;
end
另一个文件:文件名是f.m
function z=f(x,y)
z=(-y+x^2+4*x-1)/2;
end
提问者：小三3236-四级
最佳答案
你那行
y(1)=y1直接改成初值。
即x=a的时候，y的值
然后将你的for j=2:n
改成for i=2:n
应当就可以了。
哦，这个不用程序的，在MATLAB里直接调用0DE45程序就可以求解的：[t,y]=ode45(@方程[t0 t1] [初值])或者你用help ode45命令查询用法

四阶龙格库塔法求解微分方程组

四阶龙格库塔法求解微分方程组
四阶龙格库塔法是一种数值解微分方程组的方法。

它的基本思想
是将微分方程组中的每个方程都离散化，并使用预测-校正的手段来逐
步计算方程组的数值解。

具体地，四阶龙格库塔法可以分为以下步骤：
1. 对微分方程组进行离散化，即将其转化为一组差分方程。

这
一步骤的详细方法因具体的微分方程组而异，需要根据情况进行逐步
求解。

2. 使用四个预测基于中间值的校正步骤，来计算方程组的数值解。

具体地，假设当前时刻为t，并已知方程组在时刻t的解为y(t)，则可以按照以下步骤逐步计算y(t+Δt)：
（1）预测y(t+Δt)的初值为y1 = y(t) + (Δt/2)k1
（2）计算y1对应的导数k2 = f(t+Δt/2, y1)
（3）校正y(t+Δt)的初值为y2 = y(t) + (Δt/2)k2
（4）计算y2对应的导数k3 = f(t+Δt/2, y2)
（5）再次校正y(t+Δt)的初值为y3 = y(t) + Δt k3
（6）计算y3对应的导数k4 = f(t+Δt, y3)
（7）利用k1, k2, k3, k4来计算y(t+Δt)，即y(t+Δt) = y(t) + (Δt/6)(k1+2k2+2k3+k4)
3. 重复第二步，直到计算出方程组在终止时刻t1的解y(t1)。

重要的是控制步长，以避免数值计算精度不足或者计算速度过慢。

综上所述，四阶龙格库塔法是一种比较常用的解微分方程组的数
值算法。

但由于龙格库塔法的计算步骤较多，且需要进行一定的离散
化处理，因此具有一定的计算复杂度和技巧难度。

四阶龙格库塔法求解微分方程组

四阶龙格库塔法求解微分方程组微分方程是数学中的一个重要分支，它描述了自然界中许多现象的发展规律。

在实际应用中，我们经常需要求解微分方程组，以得到系统的演化规律和性质。

本文将介绍一种常用的求解微分方程组的数值方法——四阶龙格库塔法。

一、微分方程组的数值解法微分方程组是形如下式的方程集合：begin{cases} frac{dy_1}{dx}=f_1(x,y_1,y_2,cdots,y_n) frac{dy_2}{dx}=f_2(x,y_1,y_2,cdots,y_n) cdotsfrac{dy_n}{dx}=f_n(x,y_1,y_2,cdots,y_n) end{cases} 其中，$y_1,y_2,cdots,y_n$是关于自变量$x$的未知函数，$f_1,f_2,cdots,f_n$是已知的函数。

求解微分方程组的数值方法主要有以下两种：1.欧拉法欧拉法是最简单的数值方法之一，它的基本思想是将微分方程组中的导数用差分代替，从而得到一个递推公式。

具体而言，欧拉法的递推公式为：$$y_{i+1}=y_i+hf(x_i,y_i)$$其中，$h$是步长，$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。

欧拉法的精度较低，容易出现数值误差，但是它的计算速度快，适用于一些简单的微分方程组。

2.龙格库塔法龙格库塔法是一种常用的高精度数值方法，它的基本思想是将微分方程组中的导数用一系列加权的差分代替，从而得到一个递推公式。

其中，四阶龙格库塔法是最为常用的一种。

具体而言，四阶龙格库塔法的递推公式为：begin{aligned} k_1&=hf(x_i,y_i)k_2&=hf(x_i+frac{h}{2},y_i+frac{k_1}{2})k_3&=hf(x_i+frac{h}{2},y_i+frac{k_2}{2})k_4&=hf(x_i+h,y_i+k_3)y_{i+1}&=y_i+frac{k_1+2k_2+2k_3+k_4}{6} end{aligned} 其中，$k_1,k_2,k_3,k_4$是加权的差分，$h$是步长，$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。

四阶Runge-Kutta方法

实验题目3 四阶Runge-Kutta 方法摘要一阶常微分方程初值问题（6.1）0(,)()dyf x y dxy x y⎧=⎪⎨⎪=⎩的数值解法是近似计算中很重要的部分。

常微分方程初值问题的数值解法是求方程（6.1）的解在点列1(0,1,)n n n x x h n -=+= 上的近似值，这里是到的步长，一般略去下标记为。

n y n h 1n x -n x h常微分方程初值问题的数值解法一般分为两大类：（1）单步法：这类方法在计算时，只用到、和，即前一步的值。

因此，在n y 1n x +n x n y 有了初值以后就可以逐步往下计算。

典型方法如龙格–库塔方法。

()R K - （2）多步法：这类方法在计算时，除用到、和以外，还要用1n y +1n x +n x n y ，即前面步的值。

典型方法如Adams 方法。

(1,2,,;0)n p y p k k -=> k经典的方法是一个四阶的方法，它的计算公式是：R K - （6.2）112341213243(22)6(,)(,)22(,)22(,)n n n n n nn 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 +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前R K -进一步需要计算四次函数值。

f 前言利用四阶龙格-库塔方法求解微分方程的初值问题程序设计流程龙格-库塔法流程图问题1（1）TestRK4('ode1', 1, [0 -1], 5, inline('-x-1'))TestRK4('ode1', 1, [0 -1], 10, inline('-x-1'))TestRK4('ode1', 1, [0 -1], 20, inline('-x-1'))（2）TestRK4('ode2', 1, [0 1], 5, inline('1./(x+1)'))TestRK4('ode2', 1, [0 1], 10, inline('1./(x+1)'))TestRK4('ode2', 1, [0 1], 20, inline('1./(x+1)'))问题2（1）TestRK4('ode3', 3, [1 0], 5, inline('x.^2.*(exp(x)-x)'))TestRK4('ode3', 3, [1 0], 10, inline('x.^2.*(exp(x)-x)'))TestRK4('ode3', 3, [1 0], 20, inline('x.^2.*(exp(x)-x)'))（2）TestRK4('ode4', 3, [1 -2], 5, inline('2*x./(1-2*x)'))TestRK4('ode4', 3, [1 -2], 10, inline('2*x./(1-2*x)'))TestRK4('ode4', 3, [1 -2], 20, inline('2*x./(1-2*x)'))问题3（1）TestRK4('ode5', 1, [0 1/3], 5, inline('x.^2+1/3*exp(-20*x)'))TestRK4('ode5', 1, [0 1/3], 10, inline('x.^2+1/3*exp(-20*x)'))TestRK4('ode5', 1, [0 1/3], 20, inline('x.^2+1/3*exp(-20*x)'))（2）TestRK4('ode6', 1, [0 1], 5, inline('exp(-20*x)+sin(x)'))TestRK4('ode6', 1, [0 1], 10, inline('exp(-20*x)+sin(x)'))TestRK4('ode6', 1, [0 1], 20, inline('exp(-20*x)+sin(x)'))（3）TestRK4('ode7', 1, [0 0], 5, inline('exp(x).*sin(x)'))TestRK4('ode7', 1, [0 0], 10, inline('exp(x).*sin(x)'))TestRK4('ode7', 1, [0 0], 20, inline('exp(x).*sin(x)'))实验所用函数function [x,y] = RK4ODE(fun, xEnd, ini, h)% RK4ODE 用四阶Runge-Kutta法解初值问题dy/dx = f(x,y),y(x0) = y0，在x处y的值%% Synopsis: [x,y] = RK4ODE(fun, xEnd)% [x,y] = RK4ODE(fun, xEnd, ini)% [x,y] = RK4ODE(fun, xEnd, ini, h)%% Input: fun = (string) 初值问题的函数% xEnd = 使用Euler法的截止点% ini = (optional)初始条件[x0 y0]，默认为[0 0]% h = (optional)步长，默认为0.05%% Output: y = 初值问题在x处y的近似值if nargin < 3ini = [0 0]; %若未给初始条件，将初始条件设为[0 0]endif nargin < 4h = 0.05; %若未给出步长，将步长设为0.05endini = ini(:); %将初始条件转为列向量，便于判断是否正确[m,n] = size(ini);if m ~= 2 | n~= 1error('初始值必须是一个含两个元素的向量[x0 y0]');endx0 = ini(1); %初始化x0y0 = ini(2); %初始化y0x = (x0:h:xEnd)'; %构建x向量y = y0*ones(length(x), 1); %初始化y向量for j=2:length(x)k1 = h * feval(fun, x(j-1), y(j-1)); %三阶Runge-Kutta法的递推公式：y(n+1) = y(n) + (k1 + 2*k2 + 2*k3 + k4) / 6k2 = h * feval(fun, x(j-1)+h/2, y(j-1)+k1/2); % k1= h * f( x(n), y(n) )k3 = h * feval(fun, x(j-1)+h/2, y(j-1)+k2/2); % k2= h * f( x(n)+h/2, y(n)+k1/2 )k4 = h * feval(fun, x(j-1)+h, y(j-1)+k3); % k3= h * f( x(n)+h/2, y(n)+k2/2 )y(j) = y(j-1) + (k1+2*k2+2*k3+k4)/6; % k4 = h * f( x(n)+h, y(n)+k3 )endfunction TestRK4(fun, xEnd, ini, N, result)h = (xEnd - ini(1))/N;[x,y] = RK4ODE(fun, xEnd, ini, h);y0 = feval(result, x);plot(x,y,'ro',x,y0,'b-');function dydx = ode1(x,y)dydx = x+y;function dydx = ode2(x,y)dydx = -y.^2;function dydx = ode3(x,y)dydx = 2*y./x + x.^2.*exp(x);function dydx = ode4(x,y)dydx = (y.^2 + y)./x;function dydx = ode5(x,y)dydx = -20*(y-x.^2) + 2*x;function dydx = ode6(x,y)dydx = -20*y + 20*sin(x) + cos(x);function dydx = ode7(x,y)dydx = -20*(y - exp(x).*sin(x)) + exp(x).*(sin(x) + cos(x));思考题1、实验一中(1)数值解和解析解相同，(2)数值解和解析解稍有不同，因为四阶Runge-Kutta方法是以小段的线性算法来近似获得微分方程的数值解，(1)的准确解是1阶的，(2)的准确解是无限阶的，因此对于(1)数值解和解析解相同。

经典四阶龙格库塔法解一阶微分方程组

数值计算课程设计1. 经典四阶龙格库塔法解一阶微分方程组1.1 运用四阶龙格库塔法解一阶微分方程组算法分析x k 1 x k h(f 1 2 f 2 2 f 3 f 4)6hy k 1 y k h 6(g 1 2g 2 2g 3 g 4 )t k 1 t k h经过循环计算由 t 0,x 0, y 0推得 t 1,x 1,y1 t 2,x 2,y 2⋯⋯每个龙格 -库塔方法都是由一个合适的泰勒方法推导而来，使得其最终全局误差为 O h N, 一种折中方法是每次进行若干次函数求值，从而省去高阶导数计算。

4阶龙格-库塔方法 (RK4)是最常用的，它适用于一般的应用，因为它非常精准，稳定，且易于编程。

f 1 f (t k , x k , y k ) ，hf 2 f (t k h 2,x khh2f 1,y k 2g 1)hf 3 f (t k h 2,x khf 2,y k 2g 2)f 4 f (t k h ,x k hf 3,y k hg 3)g 1 g(t k ,x k , y k ) h g 2g(t k 2,x khh2f 1,y k 2g 1)hg 3 g(t k 2, x khf 2h2,y k 2g 2) g 4 g(t k h,x k hf 3, y k hg 3)1-1)1-2)1-3)1-4)1-5)1-6)1-7)1-8)1-9)1-10 )经典四阶龙格库塔法解一阶微分方程组1.2 经典四阶龙格库塔法解一阶微分方程流程图1.3 经典四阶龙格库塔法解一阶微分方程程序代码：#include <iostream> #include <iomanip> using namespace std;void RK4( double (*f)(double t,double x, double y),double (*g)(double t,double x, double y) ,double initial[3], double resu[3],double h){double f1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1; t0=initial[0];x0=initial[1];y0=initial[2]; f1=f(t0,x0,y0); g1=g(t0,x0,y0); f2=f(t0+h/2, x0+h*f1/2,y0+h*g1/2);g2=g(t0+h/2,x0+h*f1/2,y0+h*g1/2); f3=f(t0+h/2,x0+h*f2/2,y0+h*g2/2);g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2);f4=f(t0+h, x0+h*f3,y0+h*g3); g4=g(t0+h, x0+h*f3,y0+h*g3);图 1-1 经典四阶龙格库塔法解一阶微分方程流程图数值计算课程设计x1=x0+h*(f1+2*f2+2*f3+f4)/6; y1=y0+h*(g1+2*g2+2*g3+g4)/6;resu[0]=t0+h;resu[1]=x1;resu[2]=y1;}int main(){double f(double t,double x, double y);double g(double t,double x, double y);double initial[3],resu[3];double a,b,H;double t,step;int i;cout<<" 输入所求微分方程组的初值t0,x0,y0:"; cin>>initial[0]>>initial[1]>>initial[2]; cout<<" 输入所求微分方程组的微分区间[a,b]:"; cin>>a>>b;cout<<" 输入所求微分方程组所分解子区间的个数step:"; cin>>step;cout<<setiosflags(ios::right)<<setiosflags(ios::fixed)<<setprecision( 10);H=(b-a)/step;cout<< initial[0]<<setw(18)<<initial[1]<<setw(18)<<initial[2]<<endl;for(i=0;i<step;i++){ RK4( f,g ,initial, resu,H);cout<<resu[0]<<setw(20)<<resu[1]<<setw(20)<<resu[2]<<endl;initial[0]=resu[0];initial[1]=resu[1];initial[2]=resu[2];}return(0);}double f(double t,double x, double y){double dx;经典四阶龙格库塔法解一阶微分方程组dx=x+2*y; return(dx);}double g(double t,double x, double y) {double dy;dy=3*x+2*y; return(dy);}1.4 经典四阶龙格库塔法解一阶微分方程程序调试结果图示：图 1-2 经典四阶龙格库塔法解一阶微分方程算法程序调试图2. 高斯列主元法解线性方程组2.1 高斯列主元法解线性方程组算法分析使用伪代码编写高斯消元过程：for k=1 to n-1 dofor i=k+1 to n l<=a(i,k)/a(k,k) for j=k to n do a(i,j)<=a(i,j)-l*a(k,j) end %end of for j b(i)<=b(i)-l*b(k) end %end offor i end %end of for k最后得到 A ，b 可以构成上三角线性方程组接着使用回代法求解上三角线性方程组应用所编写程序计算所给例题：其中初值为求解区间为。

四阶龙格-库塔法求解常微分方程的初值问题

实验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。

四阶龙格库塔法微分方程求解发散

四阶龙格库塔法微分方程求解发散下载提示：该文档是本店铺精心编制而成的，希望大家下载后，能够帮助大家解决实际问题。

文档下载后可定制修改，请根据实际需要进行调整和使用，谢谢!本店铺为大家提供各种类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，想了解不同资料格式和写法，敬请关注!Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you! In addition, this shop provides you with various types of practical materials, such as educational essays, diary appreciation, sentence excerpts, ancient poems, classic articles, topic composition, work summary, word parsing, copy excerpts, other materials and so on, want to know different data formats and writing methods, please pay attention!四阶龙格库塔法微分方程求解发散引言微分方程是数学中一种重要的工具，用于描述自然界中许多现象的变化规律。

数值计算：四阶龙格-库塔法for二阶微分方程

数值计算：四阶龙格-库塔法for⼆阶微分⽅程引⾔考虑存在以下⼆阶偏微分⽅程f2⋅¨X(t)+f1⋅˙X(t)+f0⋅X(t)=F(t)如何使⽤四阶龙格-库塔法求解该微分⽅程？⼀阶微分⽅程的解法⾸先回顾下对于⼀阶微分⽅程的解法，现在有以下⼀阶常系数⾮齐次微分⽅程f1⋅˙X(t)+f0⋅X(t)=F(t)⽅程可以写作˙X(t)=F(t)−f0⋅X(t)f1X n+1−X ndt=F(t)−f0⋅X(t)f1X n+1=X n+dt⋅F(t)−f0⋅X nf1=X n+dt⋅f(t,X n)其中dt为时间步长。

四阶龙格-库塔法公式如下：X n+1=X n+dt6k1+2k2+2k3+k4k1=f t,X nk2=f t+dt2,Xn+dt2⋅k1k3=f t+dt2,Xn+dt2⋅k2k4=f t+dt,X n+dt⋅k3利⽤以上公式即可显式推进求解⼆阶微分⽅程的解法现在考虑⼆阶常系数⾮齐次微分⽅程f2⋅¨X(t)+f1⋅˙X(t)+f0⋅X(t)=F(t)⽅程中出现⼆次导数项¨X(t)，难以处理，所以考虑将⼆次导数项¨X(t)转化为⼀次导数项，⾄此，我们引⼊a(t)=X(t)b(t)=˙a(t)=˙X(t)原⽅程可以写成⼀阶偏微分⽅程组{()()()()(){f 2⋅˙b (t )+f 1⋅b (t )+f 0⋅a (t )=F (t )b (t )=˙a(t )仿照以上⼀阶微分⽅程的RK4解法˙a(t )=b (t )˙b (t )=F (t )−(f 1⋅b (t )+f 0⋅a (t ))f 2˙a n +1=a n +dt ⋅b (t )=a n +dt ⋅f (t ,a n ,b n )˙bn +1=b n +dt ⋅F (t )−(f 1⋅b (t )+f 0⋅a (t ))f 2=b n +dt ⋅g (t ,a n ,b n )˙a n +1=a n +dt ⋅f (t ,a n ,b n )˙bn +1=b n +dt ⋅g (t ,a n ,b n )化简⾄此，便可以仿照之前的RK4⽅法进⾏求解，具体公式为a n +1b n +1=a n +dt6⋅(k 1a +2k 2a +2k 3a +k 4a )b n +dt6⋅(k 1b +2k 2b +2k 3b +k 4b )其中，各个系数求解公式为：k 1a k 1b=f (t ,a n ,b n )g (t ,a n ,b n )k 2a k 2b =f (t +dt2,a n +dt2⋅k 1a ,b n +dt2⋅k 1b )g (t +dt2,a n +dt2⋅k 1a ,b n +dt2⋅k 1b )k 3a k 3b=f (t +dt2,a n +dt2⋅k 2a ,b n +dt2⋅k 2b )g (t +dt2,a n +dt2⋅k 2a ,b n +dt2⋅k 2b )k 4a k 4b=f (t +dt ,a n +dt ⋅k 3a ,b n +dt ⋅k 3b )g (t +dt ,a n +dt ⋅k 3a ,b n +dt ⋅k 3b )⾄此已经完成使⽤四阶龙格-库塔法⼆阶微分⽅程求解过程的梳理实例求解，以下偏微分⽅程：\begin{align} 1 \cdot \ddot{X(t)}+0 \cdot \dot{X(t)} +0 \cdot {X(t)} =cos(t)\\ \dot{X(0)}=0\\ {X(0)} =0 \end{align}理论解：\begin{align} {X(t)} =-cos(t)+1 \end{align}#include"RK_ODE.h"#include<iostream>using namespace std;#define M 1{{{{[][]{[][][][][][][][]MatrixXd F2(double t) {MatrixXd res = MatrixXd::Identity(M, M);return res;}MatrixXd F1(double t) {MatrixXd res = MatrixXd::Zero(M, M);return res;}MatrixXd F0(double t) {MatrixXd res = MatrixXd::Zero(M, M);return res;}MatrixXd F(double t) {MatrixXd res = MatrixXd::Ones(M, 1);res(0, 0) = cos(t);return res;}int main() {RK_ODE ode(M, 0.1);ode.X.fill(0.);ode.dX.fill(0.);for (int i = 0; i < 1000; i++) {ode.Solve(F2, F1, F0, F);ode.WriteMotionAndForce("results1.txt", 0); }//ode.SaveSnapshoot("snapshoot.json");cout << ode._Mode << endl;}Processing math: 88%。

matlab用四阶龙格库塔函数求解微分方程组

一、介绍Matlab作为一种强大的科学计算软件，提供了众多函数和工具来解决微分方程组。

其中，四阶龙格库塔函数是一种常用的数值方法，用于求解常微分方程组。

本文将介绍如何使用Matlab中的四阶龙格库塔函数来求解微分方程组，并对该方法的原理和实现进行详细说明。

二、四阶龙格库塔方法四阶龙格库塔方法是一种常用的数值方法，用于求解常微分方程组。

它是一种显式的Runge-Kutta方法，通过逐步逼近微分方程的解，在每一步使用多个中间值来计算下一步的解。

该方法通过四个中间值来计算下一步的状态，并且具有较高的精度和稳定性。

三、在Matlab中使用四阶龙格库塔方法求解微分方程组在Matlab中，可以使用ode45函数来调用四阶龙格库塔方法来解决微分方程组的问题。

ode45函数是Matlab提供的用于求解常微分方程组的函数，可以通过指定微分方程组以及初值条件来调用四阶龙格库塔方法来进行求解。

1. 定义微分方程组我们需要定义要求解的微分方程组。

可以使用Matlab中的匿名函数来定义微分方程组，例如：```matlabf = (t, y) [y(2); -sin(y(1))];```其中，f是一个匿名函数，用于表示微分方程组。

在这个例子中，微分方程组是y' = y2, y2' = -sin(y1)。

2. 指定初值条件和求解区间接下来，我们需要指定微分方程组的初值条件和求解区间。

初值条件可以通过指定一个初始时刻的状态向量来完成，例如：```matlabtspan = [0, 10];y0 = [0, 1];```其中，tspan表示求解区间，y0表示初值条件。

3. 调用ode45函数进行求解我们可以通过调用ode45函数来求解微分方程组的数值解。

具体的调用方式如下：```matlab[t, y] = ode45(f, tspan, y0);```其中，t和y分别表示求解的时间点和对应的状态值。

四、示例下面我们通过一个具体的例子来演示如何使用Matlab中的四阶龙格库塔方法来求解微分方程组。

《微分方程的数值解法maab四阶龙格—库塔法》PPT模板课件

k1 f (tn , yn )
k2
f (tn
1 2
h,
yn
h 2 k1)
k3
f (tn
1 2
h,
yn
h 2
k2
)
k 4 f (tn h , y n hk 3 )
四阶 Runge-Kutta 法计算流程图
开始
h 初始条件：t
迭代次数：
；y
0
N
0
积分步长：
tn t0
for i = 1 : N
解析解： x x x1 3 2(((ttt))) 0 .0 8 1 1 2 P 8 k 0siw n t) (2 .6 3 0 3 3 P k 0siw n t) (0 .2 12 2 2 P k 0siw n t)(
第一个质量的位移响应时程
各种solver 解算指令的特点
解法指令解题类型
特点
ode45 非刚性采用4、5阶Runge－Kutta法
适合场合大多数场合的首选算法
ode23 非刚性采用Adams算法
较低精度（10－3）场合
ode113 非刚性
ode23t ode15s
适度刚性
刚性
ode23s 刚性 ode23tb 刚性
y2
(0)
(2.2) (2.3)
例：著名的Van der Pol方程
y (y21)y y0
令
y1y,y2y
Y
y y
1 2
降为一阶 Yyy12(y12y12)y2y1
初始条件
Y0
y1(0) y2(0)
y10 y20
3. 根据式(2.2)编写计算导数的M函数文件ODE文件

四阶龙格库塔法(Runge-Kutta)求解微分方程

四阶龙格库塔法（Runge-Kutta ）求解微分方程张晓颖（天津大学材料学院学号：1012208027）1 引言计算传热学中通常需要求解常微分方程。

这类问题的简单形式如下：{),(')(00y x f y y x y == （1）虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中的多数微分方程需要采用数值解法求解。

初值问题（1）的数值解法有个基本特点，它们采取“步进式”，即求解过程顺着节点排序一步一步向前推进。

这类算法是要给出用已知信息y n 、 y n −i ……计算y n +1的递推公式即可。

2 龙格库塔法（Runge-Kutta ）介绍假设对于初值问题（1）有解 y = y (x ) ，用 Taylor 展开有：......)(!3)(!2)()()(321+'''+''+'+=+n n n n n x y h x y h x y h x y x y （2）龙格库塔法（Runge-Kutta ）实质上是间接的使用 Taylor 级数法的一种方法。

对于差商hx y x y n n )()(1-+，根据微分中值定理，存在 0 < θ < 1 ，使得：)()()(1h x y hx y x y n n θ+'=-+ （3）于是对于 y = y (x ) ，可得到：))(,()()(1h x y h x hf x y x y n n n n θθ+++=+ （4）设))(,(*h x y h x f K n n θθ++=，为区间 [x n , x n +1 ] 上的平均斜率。

四阶龙格库塔格式中的*K 由下式计算得到：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()2,2()2,2(),()22(6342312143211hK y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n nn n n n n （5）四阶龙格库塔法（Runge-Kutta ）的每一步需要四次计算函数值f ，其截断误差更低，计算的精度更高。