微分方程的数值解法matlab(四阶龙格—库塔法).

合集下载

matlab四阶龙格库塔法解方程组

matlab四阶龙格库塔法解方程组摘要：1.MATLAB 与四阶龙格- 库塔法简介2.四阶龙格- 库塔法求解方程组的原理3.MATLAB 中实现四阶龙格- 库塔法的方法4.四阶龙格- 库塔法在MATLAB 中的应用实例5.总结与展望正文：一、MATLAB 与四阶龙格- 库塔法简介MATLAB 是一种广泛应用于科学计算、数据分析和可视化的编程语言，它为用户提供了丰富的函数库和工具箱，使得各种数学运算和工程计算变得简单易行。

在MATLAB 中，求解方程组是工程和技术领域中常见的问题，而四阶龙格- 库塔法（RK4）是一种高效的数值求解方法。

二、四阶龙格- 库塔法求解方程组的原理四阶龙格- 库塔法是一种基于分步法的四阶数值积分方法，用于求解常微分方程初值问题。

它通过将求解区间分为若干个小区间，然后在每个小区间内，对导数进行四次评估，最后以加权平均的方式获取区间内函数的平均斜率，从而近似求得该区间内函数的值。

通过这种方式，可以逐步求解出方程组的解。

三、MATLAB 中实现四阶龙格- 库塔法的方法在MATLAB 中，可以使用自定义函数和循环结构实现四阶龙格- 库塔法求解方程组。

以下是一个简单的示例：```matlabfunction dXdt = rk4(t, X, f, dt)% 计算k1k1 = f(t, X);% 计算k2k2 = f(t + dt/2, X + 0.5*dt*k1);% 计算k3k3 = f(t + dt/2, X + 0.5*dt*k2);% 计算k4k4 = f(t + dt, X + dt*k3);% 计算四阶龙格- 库塔法导数dXdt = (k1 + 2*k2 + 2*k3 + k4) / 6;end```四、四阶龙格- 库塔法在MATLAB 中的应用实例假设我们要求解如下方程组：```x" = 2*y - zy" = x + 2*zz" = -x + y```我们可以使用MATLAB 中的四阶龙格- 库塔法求解该方程组，具体步骤如下：1.定义方程组的函数形式：```matlabfunction f = example_function(t, X)f(1, X) = [2*X(2) - X(3); X(1) + 2*X(3); -X(1) + X(2)];end```2.设置求解参数：```matlabtspan = [0, 10];dt = 0.01;```3.初始化解：```matlabX0 = [1; 1; 1];```4.使用四阶龙格- 库塔法求解方程组：```matlab[~, X] = ode45(@(t, X) example_function(t, X), tspan, X0, dt);```5.绘制解的曲线：```matlabplot3(X(:, 1), X(:, 2), X(:, 3));xlabel("x");ylabel("y");zlabel("z");title("四阶龙格- 库塔法求解示例");```五、总结与展望四阶龙格- 库塔法作为一种高效的数值积分方法，在MATLAB 中得到了广泛的应用。

matlab龙格库塔法程序,给出实例

一、介绍龙格库塔法龙格库塔法（Runge-Kutta method）是一种数值计算方法，用于求解常微分方程的数值解。

它通过多步迭代的方式逼近微分方程的解，并且具有较高的精度和稳定性。

二、龙格库塔法的原理龙格库塔法采用迭代的方式来逼近微分方程的解。

在每一步迭代中，计算出当前时刻的斜率，然后根据这个斜率来求解下一个时刻的值。

通过多步迭代，可以得到微分方程的数值解。

三、龙格库塔法的公式龙格库塔法可以表示为以下形式：k1 = f(tn, yn)k2 = f(tn + h/2, yn + h/2 * k1)k3 = f(tn + h/2, yn + h/2 * k2)k4 = f(tn + h, yn + h * k3)yn+1 = yn + h/6 * (k1 + 2k2 + 2k3 + k4)其中，k1、k2、k3、k4为斜率，h为步长，tn为当前时刻，yn为当前时刻的解，yn+1为下一个时刻的解。

四、使用matlab实现龙格库塔法在MATLAB中，可以通过编写函数来实现龙格库塔法。

下面是一个用MATLAB实现龙格库塔法的简单例子：```matlabfunction [t, y] = runge_kutta(f, tspan, y0, h)t0 = tspan(1);tf = tspan(2);t = t0:h:tf;n = length(t);y = zeros(1, n);y(1) = y0;for i = 1:n-1k1 = f(t(i), y(i));k2 = f(t(i) + h/2, y(i) + h/2 * k1);k3 = f(t(i) + h/2, y(i) + h/2 * k2);k4 = f(t(i) + h, y(i) + h * k3);y(i+1) = y(i) + h/6 * (k1 + 2*k2 + 2*k3 + k4);endend```以上就是一个简单的MATLAB函数，可以利用该函数求解给定的微分方程。

matlab代码实现四阶龙格库塔求解三阶微分方程

【Matlab代码实现四阶龙格库塔求解三阶微分方程】1. 引言在数学和工程领域中，微分方程是一种常见的数学工具，用于描述自然和工程系统的动力学行为。

而龙格-库塔（Runge-Kutta）方法是一种数值求解微分方程的常见方法。

在本文中，我将使用Matlab来实现四阶龙格-库塔方法，用以求解一个三阶微分方程的数值解。

2. 四阶龙格-库塔方法简介四阶龙格-库塔方法是一种常见的数值积分方法，用于求解常微分方程的初值问题。

该方法通过逐步逼近真实解，以得到一系列的数值解。

在Matlab中，我们可以使用内置的ode45函数来实现四阶龙格-库塔方法。

然而，在本文中，我将展示如何手动实现该方法，并解决一个具体的三阶微分方程。

3. 三阶微分方程的建立假设我们要解决的三阶微分方程为：\[y''' = f(x, y, y', y'')\]其中，\(f(x, y, y', y'')\)是一个关于自变量x和因变量y及其导数的函数。

我们的目标是使用四阶龙格-库塔方法，求解该微分方程的数值解。

4. Matlab代码实现我们需要将三阶微分方程转化为一组一阶微分方程。

设:\[z_1 = y\]\[z_2 = y'\]\[z_3 = y''\]这样，原始的三阶微分方程可以被转化为如下的一组一阶微分方程：\[z'_1 = z_2\]\[z'_2 = z_3\]\[z'_3 = f(x, z_1, z_2, z_3)\]接下来，我们可以使用四阶龙格-库塔方法来求解上述一组一阶微分方程。

以下是我在Matlab中手动实现四阶龙格-库塔方法的代码：```matlabfunction [t, y] = myRK4(f, a, b, y0, h)N = (b - a) / h;t = zeros(N+1, 1);y = zeros(N+1, length(y0));t(1) = a;y(1, :) = y0;for i = 1:Nk1 = h * f(t(i), y(i, :));k2 = h * f(t(i) + h/2, y(i, :) + k1/2);k3 = h * f(t(i) + h/2, y(i, :) + k2/2);k4 = h * f(t(i) + h, y(i, :) + k3);y(i+1, :) = y(i, :) + (k1 + 2*k2 + 2*k3 + k4) / 6;t(i+1) = t(i) + h;endend```在上述代码中，我定义了一个名为myRK4的函数，用于实现四阶龙格-库塔方法。

MATLAB常微分方程数值解——欧拉法、改进的欧拉法与四阶龙格库塔方法

MATLAB常微分⽅程数值解——欧拉法、改进的欧拉法与四阶龙格库塔⽅法MATLAB常微分⽅程数值解作者：凯鲁嘎吉 - 博客园1.⼀阶常微分⽅程初值问题2.欧拉法3.改进的欧拉法4.四阶龙格库塔⽅法5.例题⽤欧拉法，改进的欧拉法及4阶经典Runge-Kutta⽅法在不同步长下计算初值问题。

步长分别为0.2,0.4,1.0.matlab程序：function z=f(x,y)z=-y*(1+x*y);function R_K(h)%欧拉法y=1;fprintf('欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K=f(x,y);y=y+h*K;fprintf('欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%改进的欧拉法y=1;fprintf('改进的欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h,y+h*K1);y=y+(h/2)*(K1+K2);fprintf('改进的欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%龙格库塔⽅法y=1;fprintf('龙格库塔法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h/2,y+(h/2)*K1);K3=f(x+h/2,y+(h/2)*K2);K4=f(x+h,y+h*K3);y=y+(h/6)*(K1+2*K2+2*K3+K4);fprintf('龙格库塔法：x=%f, y=%f\n',x+h,y);end结果：>> R_K(0.2)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.200000, y=0.800000欧拉法：x=0.400000, y=0.614400欧拉法：x=0.600000, y=0.461321欧拉法：x=0.800000, y=0.343519欧拉法：x=1.000000, y=0.255934改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.200000, y=0.807200改进的欧拉法：x=0.400000, y=0.636118改进的欧拉法：x=0.600000, y=0.495044改进的欧拉法：x=0.800000, y=0.383419改进的欧拉法：x=1.000000, y=0.296974龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.200000, y=0.804636龙格库塔法：x=0.400000, y=0.631465龙格库塔法：x=0.600000, y=0.489198龙格库塔法：x=0.800000, y=0.377225龙格库塔法：x=1.000000, y=0.291009>> R_K(0.4)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.400000, y=0.600000欧拉法：x=0.800000, y=0.302400改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.400000, y=0.651200改进的欧拉法：x=0.800000, y=0.405782龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.400000, y=0.631625龙格库塔法：x=0.800000, y=0.377556>> R_K(1)欧拉法：x=0.000000, y=1.000000欧拉法：x=1.000000, y=0.000000改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=1.000000, y=0.500000龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=1.000000, y=0.303395注意：在步长h为0.4时，要将for i=1:1/h改为for i=1:0.8/h。

微分方程的数值解法matlab(四阶龙格—库塔法).ppt

3月15日作业: 1.Van der Pol 方程的两种解法:1)采用ode45命令 2)Runge-Kutta方法 2.Duffing 方程的求解（Runge-Kutta方法，计算步长h=0.005，计
算时间t0=0.0,tN=100）
要求：写出程序体，打印所绘图形，图形标题用个人的名字。
Duffing 方程Fra bibliotek子程序：RK_sub.m function ydot = vdpol (t, y) ydot=zeros(size(y)); ydot(1) = y(2); ydot(2) = -y(2)*(y(1)^2-1)-y(1); 或写为：
ydot = [y(1) ;-y(2)*(y(1)^2-1)-y(1)];
K4 = ZCX_sub (t0 + h, y0 + hK3, P) Y1 = y0 + (h/6) (K1 + 2K2 + 2K3 + K4) t1, Y1 (输出 t1, y1) next i 输出数据或图形
(t ) AY (t ) P (t ) Y
Matlab 程序(子程序：ZCX_sub.m)
多步法-Admas方法
计算
的近似值 y(tn1 )
时只用到 n1
y
，是自开始方法 n n
t ,y

Runge-Kutta法是常微分方程的一种经典解法 MATLAB 对应命令：ode45
四阶Runge-Kutta公式
h yn 1 yn (k1 2k 2 2k3 k 4 ) 6 k1 f (t n , yn ) 1 h k 2 f (t n h, yn k1 ) 2 2 1 h k3 f (t n h, yn k 2 ) 2 2 k 4 f (t n h, yn hk3 )

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

函数功能编辑本段回目录ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。

ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)³。

解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解.使用方法编辑本段回目录[T,Y] = ode45(odefun,tspan,y0)odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名tspan 是区间[t0 tf] 或者一系列散点[t0,t1,...,tf]y0 是初始值向量T 返回列向量的时间点Y 返回对应T的求解列向量[T,Y] = ode45(odefun,tspan,y0,options)options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等[T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options)在设置了事件参数后的对应输出TE 事件发生时间YE 事件解决时间IE 事件消失时间sol =ode45(odefun,[t0 tf],y0...)sol 结构体输出结果应用举例编辑本段回目录1 求解一阶常微分方程程序:一阶常微分方程odefun=@(t,y) (y+3*t)/t^2; %定义函数tspan=[1 4]; %求解区间y0=-2; %初值[t,y]=ode45(odefun,tspan,y0);plot(t,y) %作图title('t^2y''=y+3t,y(1)=-2,1<t<4')legend('t^2y''=y+3t')xlabel('t')ylabel('y')% 精确解% dsolve('t^2*Dy=y+3*t','y(1)=-2')% ans =一阶求解结果图% (3*Ei(1) - 2*exp(1))/exp(1/t) - (3*Ei(1/t))/exp(1/t)2 求解高阶常微分方程关键是将高阶转为一阶,odefun的书写.F(y,y',y''...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为列向量dxdy=[y(1),y(2)....]程序:function Testode45tspan=[3.9 4.0]; %求解区间y0=[2 8]; %初值[t,x]=ode45(@odefun,tspan,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')title('y'' ''=-t*y + e^t*y'' +3sin2t')xlabel('t')ylabel('y')function y=odefun(t,x)y=zeros(2,1); % 列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);endend高阶求解结果图相关函数编辑本段回目录ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tbMatlab中龙格-库塔(Runge-Kutta)方法原理及实现(自己写的，非直接调用)龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。

matlab编写龙格库塔法或欧拉法求解常微分方程数值解

龙格库塔法（Runge-Kutta method）和欧拉法（Euler's method）是两种常用的数值求解常微分方程的方法。

这里分别给出它们的MATLAB实现：1. 龙格库塔法（Runge-Kutta method）：```matlabfunction [y, t] = runge_kutta(f, y0, t0, tf, h)% f: 微分方程函数，输入为[y, t]，输出为dy/dt% y0: 初始值% t0: 初始时间% tf: 结束时间% h: 步长N = round((tf - t0) / h); % 计算迭代次数t = zeros(1, N + 1); % 初始化时间向量y = zeros(size(y0), N + 1); % 初始化解向量t(1) = t0;y(:, 1) = y0;for i = 1:Nk1 = h * f(y(:, i), t(i));k2 = h * f(y(:, i) + k1 / 2, t(i) + h / 2);k3 = h * f(y(:, i) + k2 / 2, t(i) + h / 2);k4 = h * f(y(:, i) + k3, t(i) + h);y(:, i + 1) = y(:, i) + (k1 + 2 * k2 + 2 * k3 + k4) / 6;t(i + 1) = t(i) + h;endend```2. 欧拉法（Euler's method）：```matlabfunction [y, t] = euler_method(f, y0, t0, tf, h)% f: 微分方程函数，输入为[y, t]，输出为dy/dt% y0: 初始值% t0: 初始时间% tf: 结束时间% h: 步长N = round((tf - t0) / h); % 计算迭代次数t = zeros(1, N + 1); % 初始化时间向量y = zeros(size(y0), N + 1); % 初始化解向量t(1) = t0;y(:, 1) = y0;for i = 1:Ny(:, i + 1) = y(:, i) + h * f(y(:, i), t(i));t(i + 1) = t(i) + h;endend```使用这两个函数时，需要定义一个表示微分方程的函数`f`，例如：```matlabfunction dydt = my_ode(y, t)dydt = -y; % 一个简单的一阶线性微分方程：dy/dt = -yend```然后调用相应的求解函数，例如：```matlaby0 = 1; % 初始值t0 = 0; % 初始时间tf = 5; % 结束时间h = 0.1; % 步长[y_rk, t_rk] = runge_kutta(@my_ode, y0, t0, tf, h);[y_euler, t_euler] = euler_method(@my_ode, y0, t0, tf, h);```。

龙格-库塔法(Runge-Kutta)matlab代码及含义

龙格－库塔法（Runge-Kutta）数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。

经典四阶龙格库塔法RK4””或者就是龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4“龙格库塔法”。

令初值问题表述如下。

则，对于该问题的RK4由如下方程给出：其中这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。

该斜率是以下斜率的加权平均：k1是时间段开始时的斜率；k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn+h/2的值；k3也是中点的斜率，但是这次采用斜率k2决定y值；k4是时间段终点的斜率，其y值用k3决定。

当四个斜率取平均时，中点的斜率有更大的权值：RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。

注意上述公式对于标量或者向量函数(y可以是向量)都适用。

显式龙格库塔法显示龙格－库塔法是上述RK4法的一个推广。

它由下式给出其中(注意：上述方程在不同著述中由不同但却等价的定义)。

要给定一个特定的方法，必须提供整数s(阶段数)，以及系数aij(对于1≤j<i≤s),bi(对于i=1,2,2,...,...,s)和ci(对于i=2,3,3,...,...,s)。

这些数据通常排列在一个助记c3a31如果要求方法有精度p则还有相应的条件，也就是要求舍入误差为O(hp+1)时的条件。

这些可以从舍入误差本身的定义中导出。

例如，一个2阶精度的2段方法要求b1+b2=1,b2c2=1/2,以及b2a21=1/2。

matlab求微分方程初值问题的特解

在MATLAB中求解微分方程初值问题，一种常用的方法是使用MATLAB的内置函数ode45。

这是一个基于四阶龙格-库塔法的方法，适用于大多数初值问题。

以下是一个简单的例子来说明如何使用ode45求解微分方程初值问题。

考虑这样一个微分方程初值问题：
dy/dx = y, 当x=0时，y=1。

在MATLAB中使用ode45求解此问题的代码如下：
matlab复制代码
% 定义微分方程函数
function dy = odeExample(x,y)
dy = y;
end
% 定义初始条件
x0 = 0;
y0 = 1;
% 定义x的范围
xspan = [01];
% 使用ode45求解
[x,y] = ode45(@odeExample, xspan, y0);
% 绘制解的图形
plot(x, y)
xlabel('x')
ylabel('y')
title('Solution of dy/dx = y')
以上代码中，首先定义了一个函数odeExample，这个函数描述了微分方程的右侧，也就是dy/dx = y的部分。

然后定义了初始条件x0和y0，再通过ode45函数求解微分方程，最后使用plot函数绘制了解的图形。

请注意，对于更复杂的微分方程，你可能需要修改odeExample函数以适应不同的形式。

同时，ode45也允许你更改步长等参数以适应更复杂的问题。

如果需要了解更多关于ode45的使用，可以参考MATLAB的官方文档。

三阶、四阶龙格库塔函数matlab代码

三阶龙格-库塔法的计算公式为：)4(6)2,()2,2(),(3211213121K K K h y y hK hK y h x g K K h y h x g K y x g K i i i i i i i i +++=+-+=++==+ 三阶龙格—库塔公式的Matlab 程序代码：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;endformat short ；DELGKT3_kuta函数运行时需要调用下列函数：function fv=Funval(f, varvec, varval ）var= findsym(f ）;if length （var ）<4if var （1)==varvec （1)fv=subs （f ，varvec （1），varval(1）);elsefv=subs(f ，varvec(2)，varval （2)）;endelsefv=subs(f ，varvec ，varval ）；end三阶龙格—库塔求解一阶常微分方程应用实例。

用三阶龙格—库塔法求下面常微分方程的数值解。

⎪⎩⎪⎨⎧=+-=1)0(232y y x dx dy 10≤≤x在编辑窗口输入下列程序段,然后执行该程序.syms x y ；z=2*x-3*y+2；yy=DELGKT3_kuta(z ，0。

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中的四阶龙格库塔方法来求解微分方程组。

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

【SOS】在matlab中四阶Runge-Kutta法求解常微分方程匚悬赏分：150 -解决时间：2008-9-16 10:19【SOS】在matlab中四阶Runge-Kutta法求解常微分方程dx1/dt=x2 x1(0)=1e-8dx2/dt=x3 x2(0)=0dx3/dt=x4 x3(0)=0dx4/dt=f(t,x1,x2,x3,x4) x4(0)=0求解区间［0,1e-6］，在matlab中用四阶Runge-Kutta法求解，怎么用从一阶循环到四阶？有没有类似的具体算例？最好能帮小弟写个子程序希望各路高人出手相助，关乎小弟能否明年一月份顺利毕业，在此先多谢各位高人了。

提问者：zhihaikou -三级最佳答案没试过matlab，算这玩意太慢了，有fortran版的要不，有兴趣的话可以参考一下。

代码:SUBROUTINE run ge_kutta()!关于Runge-Kutta方法，该方法是用来解形如y'=f(t,y)的常微分方程的!经典的4阶R-K方法的公式如下：! Yn+1 = Yn + h/6 * (K1+2K2+2K3+K4)!其中! K仁f(Tn,Yn)! K2=f(T n+h/2,Y n+h/2*K1)! K3=f(T n+h/2,Y n+h/2*K2)! K4=f(T n+h,Y n+h*K3)type( no de_fish),po in ter :: pi7HO!S%（『!）siuo!」inu% 厂u 人 NOa%（r!）siuau;nu%r u 人NOd%（『!）siuo!」inu% 厂 u 人t7HN%（『!）siuo!」inu% 厂 u 人£ON%（『!）siuo!」inu%厂u 人（（>HO!S%lueu;nu _U!UJ >（ NOa%lueu;nu _u!iu > （ NOd%lueu;nu _u!iu >（ >HN%lueu;nu _U!UJ > eON%lueu;nu _u!iu >WW : >d31S -一ii（W+£>k0N+Z>k0N+"）*（09dois —|j ） + u 人二厂u 人L+u 人直具 *-以吐u 人甲：ed3丄S-—一ii（e>i*de;s -p+u 人'dois —：M+iuaurK ）—：|j ）uo!Qun 厂 |E 」60屮!=厂》£>UL|+u 人 POM u人'u+l DOMfeJB來：kSd3丄S ------- i（2>1*（0S/deis _p ）+UA 10S/deis _p+juajjno _p ）uo!pun4_|ej6aju!=e>j◊u ⑵q ）+u 人［1 變歪 u 人'乙/i|+i pOMfeJB Ux 來：c sd3丄s ------------- i（L>1*（0S/deis _|j ）+UA 107/deis _p+juajjno _p ）uo!pun4_|ej6aju!=3>j以*⑵L|）+u 人 H 變歪U 人'乙/1|+1 POMHB 辽X 來：2 2 C □丄s --------- i（uA^uejjno-pJuoipunPie 」6。

matlab用四阶龙格库塔法解二阶常微分方程

matlab用四阶龙格库塔法解二阶常微分方程在数学和工程中，常微分方程是描述自然界中各种物理现象和过程的常见数学模型。

常微分方程通常包含未知函数及其导数的方程。

在求解常微分方程时，人们通常使用数值方法来近似求解，其中一种常见的方法是使用龙格-库塔方法。

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

它是通过在每个步骤中估计未知函数的导数来数值求解方程。

它的精度比较高，通常被认为是最常用的龙格-库塔方法。

对于一个二阶常微分方程：y''(t) = f(t, y(t), y'(t))其中y(t)是未知函数，f(t, y(t), y'(t))是已知的函数。

我们可以将这个二阶微分方程转化为两个一阶微分方程：y'(t) = u(t)u'(t) = f(t, y(t), u(t))其中u(t)是y'(t)的一个替代函数。

为了使用四阶龙格-库塔方法求解这个方程，我们需要将时间区间[t0, tn]分割为等距的n个子区间，每个子区间的长度为h = (tn - t0)/n。

我们从初始点t0开始，通过多个步骤来逼近解直到tn。

每一步我们使用以下递推公式来估计y(t)和u(t)的值：k1 = h * u(t)l1 = h * f(t, y(t), u(t))k2 = h * (u(t) + l1/2)l2 = h * f(t + h/2, y(t) + k1/2, u(t) + l1/2)k3 = h * (u(t) + l2/2)l3 = h * f(t + h/2, y(t) + k2/2, u(t) + l2/2)k4 = h * (u(t) + l3)l4 = h * f(t + h, y(t) + k3, u(t) + l3)然后我们可以使用以下公式来更新y(t)和u(t)的值：y(t + h) = y(t) + (k1 + 2k2 + 2k3 + k4)/6u(t + h) = u(t) + (l1 + 2l2 + 2l3 + l4)/6这个过程重复n次，直到达到结束时间tn。

matlab算法-求解微分方程数值解和解析解

MATLAB是一种用于数学计算、工程和科学应用程序开发的高级技术计算语言和交互式环境。

它被广泛应用于各种领域，尤其在工程和科学领域中被用于解决复杂的数学问题。

微分方程是许多工程和科学问题的基本数学描述，求解微分方程的数值解和解析解是MATLAB算法的一个重要应用。

1. 求解微分方程数值解在MATLAB中，可以使用各种数值方法来求解微分方程的数值解。

其中，常见的方法包括欧拉法、改进的欧拉法、四阶龙格-库塔法等。

这些数值方法可以通过编写MATLAB脚本来实现，从而得到微分方程的近似数值解。

以常微分方程为例，可以使用ode45函数来求解微分方程的数值解。

该函数是MATLAB中用于求解常微分方程初值问题的快速、鲁棒的数值方法，可以有效地得到微分方程的数值解。

2. 求解微分方程解析解除了求解微分方程的数值解外，MATLAB还可以用于求解微分方程的解析解。

对于一些特定类型的微分方程，可以使用符号计算工具箱中的函数来求解微分方程的解析解。

通过符号计算工具箱，可以对微分方程进行符号化处理，从而得到微分方程的解析解。

这对于研究微分方程的性质和特点非常有帮助，也有助于理论分析和验证数值解的准确性。

3. MATLAB算法应用举例在实际工程和科学应用中，MATLAB算法求解微分方程问题非常常见。

在控制系统设计中，经常需要对系统的动态特性进行分析和设计，这通常涉及到微分方程的建模和求解。

通过MATLAB算法，可以对系统的微分方程进行数值求解，从而得到系统的响应曲线和动态特性。

另外，在物理学、生物学、经济学等领域的建模和仿真中，也经常需要用到MATLAB算法来求解微分方程问题。

4. MATLAB算法优势相比于其他数学软件和编程语言，MATLAB在求解微分方程问题上具有明显的优势。

MATLAB提供了丰富的数值方法和工具，能够方便地对各种微分方程进行数值求解。

MATLAB具有直观的交互式界面和强大的绘图功能，能够直观地展示微分方程的数值解和解析解，有利于分析和理解问题。

MATLAB龙格库塔法

MATLAB龙格库塔法⾮刚性常微分⽅程的数值解法通常会⽤四阶龙格库塔算法，其matlab函数对应ode45。

对于dy/dx = f(x,y)，y(0)=y0。

其四阶龙格库塔公式如下：对于通常计算，四阶已经够⽤，四阶以上函数f(x,y)计算⼯作量⼤⼤增加⽽精度提⾼较慢。

下⾯以龙格库塔法解洛伦兹⽅程为例：matlab代码如下：main.m：1 clear all;2 close all;3 clc;45 %系统龙格库塔法6 [t,h] = ode45(@test_fun,[040],[1240]);7 plot3(h(:,1),h(:,2),h(:,3));8 grid on;910 %⾃定义龙格库塔法11 [t1,h1]=runge_kutta(@test_fun,[1240],0.01,0,40);12 figure;13 plot3(h1(1,:),h1(2,:),h1(3,:),'r')14 grid on;runge_kutta.m（函数参考⽹络）：1 %参数表顺序依次是微分⽅程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）2 function [x,y]=runge_kutta(ufunc,y0,h,a,b)3 n=floor((b-a)/h); %步数4 x(1)=a; %时间起点5 y(:,1)=y0; %赋初值，可以是向量，但是要注意维数6for i=1:n %龙格库塔⽅法进⾏数值求解7 x(i+1)=x(i)+h;8 k1=ufunc(x(i),y(:,i));9 k2=ufunc(x(i)+h/2,y(:,i)+h*k1/2);10 k3=ufunc(x(i)+h/2,y(:,i)+h*k2/2);11 k4=ufunc(x(i)+h,y(:,i)+h*k3);12 y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;13 endtest_fun（洛伦兹⽅程）：1 %构造微分⽅程2 function dy=test_fun(t,y)3 a = 16;4 b = 4;5 c = 45;67 dy=[a*(y(2)-y(1));8 c*y(1)-y(1)*y(3)-y(2);9 y(1)*y(2)-b*y(3)];得到很经典的洛伦兹吸引⼦，结果如下：。

微分方程的数值解法matlab(四阶龙格—库塔法)说课材料PPT文档38页

1
0
、
倚
南
窗
以
寄
傲
，
审
容
膝
之
易
安
。
61、奢侈是舒适的，否则就不是奢侈。——CocoCha nel 62、少而好学，如日出之阳；壮而好学，如日中之光；志而好学，如炳烛之光。 ——刘向 63、三军可夺帅也，匹夫不可夺志也。 ——孔丘 64、人生就是学校。在那里，与其说好的教师是幸福，不如说好的教师是不幸。 ——海贝尔 65、接受挑战，就可以享受胜利的喜悦。——杰纳勒尔·乔治·S·巴顿
微分方程的数值解法matlab(四阶龙格—库塔法)说课材料
6
、
露
凝
无
游
氛
，
天
高
风
景
澈
。
7、翩翩新来燕，双双入我庐，先巢故尚在，相将还旧居。
ቤተ መጻሕፍቲ ባይዱ
8
、
吁
嗟
身
后
名
，
于
我
若
浮
烟
。
9、陶渊明（约 365年 —427年），字元亮，（又一说名潜，字渊明）号五柳先生，私谥“靖节”，东晋末期南朝宋初期诗人、文学家、辞赋家、散
文家。汉族，东晋浔阳柴桑人（今江西九江）。曾做过几年小官，后辞官回家，从此隐居，田园生活是陶渊明诗的主要题材，相关作品有《饮酒》、《归园田居》、《桃花源记》、《五柳先生传》、《归去来兮辞》等。
谢谢！

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

(2.2)
y1 (0) Y0 y ( 0 ) 2
(2.3)
例：著名的Van der Pol方程
2 y0 y ( y 1) y
令
y1 y, y2 y
y1 Y y2
降为一阶
1 y y 2 Y y 2 2 ( y1 1) y2 y1
四阶 Runge-Kutta 法计算流程图
开始初始条件： t 0；y0 迭代次数： N 积分步长：
h
t n t0
for i = 1 : N
y2 (t0 ) y2 (t1 ) y2 (t 2 )
Y (:,1) y1 (t ) Y (:,2) y2 (t )
例题1：著名的Van der Pol方程
2 y0 y ( y 1) y
解法1：采用 ODE命令
% 主程序（程序名：VanderPol _ex1.m） t0 = 0; tN = 20; tol = 1e-6; Y0 = [0.25; 0.0]; [t, Y]=ode45 (‘dYdt’, t0, tN, Y0, tol); subplot (121), plot (t, Y)
%M function file name: dYdt.m function Yd = f (t, Y) Yd=zeros(size(Y)); Yd (1) Y (2); Yd (2) (Y (1).^2 1) *Y (2) Y (1);
4. 使编写好的ODE函数文件和初值 Y0 供微分方程解算指令（solver）调用
ode23tb 刚性
二. 四阶 Runge-Kutta 法
dy f (t , y) a t b dt y(t0 ) y0
对 I=[a,b]作分割步长
a t0 t1 t N 1 b
hi hi 1 hi , i 0,1,, N 1
初值问题的数值解法分为两大类
y1 (0) y10 Y0 y 2 (0) y 20
初始条件
3. 根据式(2.2)编写计算导数的M函数文件ODE文件
作为输出宗量把t,Y作为输入宗量，把 Y
%M function file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展开式例Van der Pol方程
ode23 ode113
ode15s
ode23tb
ቤተ መጻሕፍቲ ባይዱ
ode23s
一.解ODE的基本机理:
1. 列出微分方程
f ( y, y , t) y
(2.1)
令
y1 y, y2 y
(0) y 0 y(0) y0 , y
初始条件
2. 把高阶方程转换成一阶微分方程组 1 y f1 (t , Y ) Y f (t , Y ) y f ( t , Y ) 2 2
说明：
Solver解算指令的使用格式
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol); ode45
输出宗量形式
y1 (t0 ) y (t ) Y 1 1 y1 (t 2 )
subplot (122), plot (Y( :, 1), Y( :, 2))
Van der Pol方程 % 子程序（程序名： dYdt.m ）
function Ydot = dYdt (t, Y) Ydot=[Y(2);-Y(2)*(Y(1)^2-1)-Y(1)];
或写为
function Ydot = dYdt (t, Y) Ydot=zeros(size(Y)); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).^2-1)-Y(1)];
单步法-Runge-Kutta 方法
多步法-Admas方法
计算 y(tn1 ) 的近似值 yn1 时只用到 t n , yn ，是自开始方法
Runge-Kutta法是常微分方程的一种经典解法
MATLAB 对应命令：ode45
四阶Runge-Kutta公式
h yn 1 yn (k1 2k 2 2k3 k 4 ) 6 k1 f (t n , yn ) 1 h k 2 f (t n h, yn k1 ) 2 2 1 h k3 f (t n h, yn k 2 ) 2 2 k 4 f (t n h, yn hk3 )
微分方程的数值解法
四阶龙格—库塔法（The Fourth-Order Runge－Kutta Method）
常微分方程（Ordinary differential equations, ODE)

初值问题---给出初始值边值问题---给出边界条件
ode45 ode23t
与初值常微分方程解算有关的指令
各种solver 解算指令的特点
解法指令 ode45 ode23 ode113 ode23t ode15s ode23s 解题类型非刚性非刚性非刚性适度刚性刚性刚性特点采用4、5阶Runge－Kutta法采用Adams算法多步法；采用Adams算法；高低精度均可（10－3～10－6）采用梯形法则算法多步法；采用2阶Rosenbrock 算式，精度中等一步法；采用2阶Rosenbrock 算式，低精度采用梯形法则－反向数值微分两阶段算法，低精度适合场合大多数场合的首选算法较低精度（10－3）场合 ode45计算时间太长时取代ode45 适度刚性当ode45失败时使用；或存在质量矩阵时低精度时，比ode15s有效；或存在质量矩阵时低精度时，比ode15s有效；或存在质量矩阵时