四阶龙格库塔法解微分方程
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]。
四阶龙格库塔实验报告
三、四阶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 。
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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
四阶龙格库塔法解微分方程
一、四阶龙格库塔法解一阶微分方程
①一阶微分方程：cos
，初始值y(0)=0，求
y t
解区间[0 10]。
MATLAB程序：
%%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程%%%%%%%%%%% y'=cost
%%%%%%%%%%% y(0)=0, 0≤t≤10，h=0.01 %%%%%%%%%%% y=sint
h=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;
end
subplot(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; end
subplot(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]。
MATLAB 程序：
%%%%%%%%% 龙格库塔欧拉方法解二阶微分方程
%%%%%%%%% y''=cost
%%%%%%%%% y'(0)=-1, y(0)=0
%%%%%%%%% 转换为一阶微分方程
h=0.01;
ht=10;
t=0:h:ht;
%%%%%%%%% y'=z1, z1'=y''=cost
y=zeros(1,length(t));
z=zeros(1,length(t));
y(1)=0;
z(1)=-1;
f=@(t,y,z)(z);
g=@(t,y,z)(sin(t));
for i=1:(length(t)-1)
K1=f(t(i),y(i),z(i));
L1=g(t(i),y(i),z(i));
K2=f(t(i)+h/2,y(i)+h/2*K1,z(i)+h/2*L1);
L2=g(t(i)+h/2,y(i)+h/2*K1,z(i)+h/2*L1);
K3=f(t(i)+h/2,y(i)+h/2*K2,z(i)+h/2*L2);
L3=g(t(i)+h/2,y(i)+h/2*K2,z(i)+h/2*L2);
K4=f(t(i)+h,y(i)+h*K3,z(i)+h*L3);
L4=g(t(i)+h,y(i)+h*K3,z(i)+h*L3);
y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4);
z(i+1)=z(i)+h/6*(L1+2*L2+2*L3+L4);
end
plot(t,y)
xlabel('t');
ylabel('y');
title('y''=sin(t)');
legend('y''=sint');
图3
②二阶微分方程：
7.35499*cos y x ，初始值y(1)=0，y'(1)=0，求解
区间[0 25]。
MATLAB程序：
%%%%%%%%% 龙格库塔欧拉方法解二阶微分方程
%%%%%%%%% y''=7.35499*cosx
%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)-1
K1=RK_z(i);
L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 and L1
K2=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;
L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4 RK_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);
end
plot(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) %周期
function w=rightf_sys1(x,y,z)
w=7.35499*cos(y);
end
图4
注：四阶龙格库塔法求解二阶微分方程时，只需把二阶微分方程转化为一阶微分方程，再进行相关求解即可。