MATLAB龙格-库塔方法解微分方程

姓名：樊元君学号：02 日期：一、实验目的掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。

掌握使用MATLAB程序求解常微分方程问题的方法。

：二、实验内容1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。

实验中以下列数据验证程序的正确性。

求，步长h=。

*2、实验注意事项的精确解为，通过调整步长，观察结果的精度的变化^)三、程序流程图：●改进欧拉格式流程图：~|●四阶龙格库塔流程图：]四、源程序：●改进后欧拉格式程序源代码：function [] = GJOL(h,x0,y0,X,Y)format longh=input('h=');…x0=input('x0=');y0=input('y0=');disp('输入的范围是：');X=input('X=');Y=input('Y=');n=round((Y-X)/h);\i=1;x1=0;yp=0;yc=0;for i=1:1:nx1=x0+h;yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);%·yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);%y1=(yp+yc)/2;x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y);:endend●四阶龙格库塔程序源代码：function [] = LGKT(h,x0,y0,X,Y)。

format longh=input('h=');x0=input('x0=');y0=input('y0=');disp('输入的范围是：');"X=input('X=');Y=input('Y=');n=round((Y-X)/h);i=1;x1=0;k1=0;k2=0;k3=0;k4=0;for i=1:1:n~x1=x0+h;k1=-x0*y0^2;%k1=y0-2*x0/y0;%k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);%…y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);%x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y);end·end*五、运行结果：改进欧拉格式结果：;}四阶龙格库塔结果：步长分别为：和时，不同结果显示验证了步长减少，对于精度的提高起到很大作用，有效数字位数明显增加。

matlab龙格库塔方法求解二元二阶常微分方程组文章标题：深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用摘要：在科学与工程领域，常常需要求解复杂的微分方程组，而matlab作为一种强大的数学工具，提供了许多求解微分方程组的方法。

本文将深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用，以便读者全面理解该方法并能灵活应用于实际问题中。

正文：一、介绍龙格库塔方法龙格-库塔法（Runge-Kutta methods）是一种数值求解常微分方程的方法，通过将微分方程的解进行离散化，将微分方程转化为差分方程，从而进行数值求解。

龙格库塔方法通过迭代计算，能够得到微分方程的数值解，广泛应用于科学计算和工程技术领域。

二、matlab中的龙格库塔方法在matlab中，龙格库塔方法通过ode45函数实现，该函数能够对一阶或高阶常微分方程进行数值求解。

用户可以通过设定初始条件、微分方程表达式，以及积分区间等参数，快速得到微分方程的数值解。

ode45函数采用自适应步长的方式进行求解，能够有效解决微分方程解的数值稳定性和精确度问题。

三、龙格库塔方法在求解二元二阶常微分方程组中的应用考虑如下形式的二元二阶常微分方程组：$$\begin{cases}y_1' = f_1(t, y_1, y_2) \\y_2' = f_2(t, y_1, y_2)\end{cases}$$其中，$y_1(t)$和$y_2(t)$是未知函数，$f_1(t, y_1, y_2)$和$f_2(t,y_1, y_2)$分别表示其对应的函数表达式。

通过matlab中的ode45函数，可以将该二元二阶常微分方程组转化为一阶常微分方程组的形式，然后利用龙格库塔方法进行数值求解。

设定初始条件$y_1(0) = y1_0, y_2(0) = y2_0$，对应的一阶方程组为：$$\begin{cases}u_1' = u_3 \\u_2' = u_4 \\u_3' = f_1(t, u_1, u_2) \\u_4' = f_2(t, u_1, u_2)\end{cases}$$其中，$u_1(t) = y_1(t), u_2(t) = y_2(t), u_3(t) = y_1'(t), u_4(t) =y_2'(t)$，通过ode45函数求解该一阶常微分方程组即可得到原二元二阶常微分方程组的数值解。

解析解： x x x1 3 2(((ttt))) 0 .0 8 1 1 2 P k 8 0siw n t) (2 .6 3 0 3 3 P k 0siw n t) (0 .2 12 2 2 P k 0siw n t)(
第一个质量的位移响应时程
Y (t)A(Y t)P(t)
(2)
Y (t)A(Y t)P(t)
3. Matlab 程序(主程序：ZCX)
t0;Y0;h;N;P0,w; %输入初始值、步长、迭代次数、初始激励力；
for i = 1 : N
t1 = t0 + h
P=[P0*sin(w*t0);0.0;0.0]
%输入t0时刻的外部激励力
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)];
Solver解算指令的使用格式
说明：
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol);
ode45
输出宗量形式
y1 (t0 )
Y
y1
(t1
)
y
1
(t
2
)
y2 (t0 )
y
2
(
t1
)
y
2
(
t

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求解常微分方程组常微分方程组是数学中的一个重要分支，它描述了多个变量随时间变化的关系。

在实际应用中，常微分方程组经常被用来描述物理、化学、生物等领域中的动态系统。

本文将介绍如何使用MATLAB求解常微分方程组。

MATLAB是一种强大的数学软件，它提供了许多工具和函数来求解常微分方程组。

在MATLAB中，我们可以使用ode45函数来求解常微分方程组。

ode45函数是一种常用的数值求解器，它使用龙格-库塔方法来求解常微分方程组。

我们需要定义常微分方程组。

常微分方程组通常采用向量形式表示，例如：dy/dt = f(t,y)其中，y是一个向量，f(t,y)是一个向量函数。

在MATLAB中，我们可以使用匿名函数来定义f(t,y)。

例如，如果我们要求解以下常微分方程组：dy1/dt = -y1 + 2*y2dy2/dt = -2*y1 + 3*y2我们可以定义f(t,y)为：f = @(t,y) [-y(1) + 2*y(2); -2*y(1) + 3*y(2)];接下来，我们需要指定初值条件。

初值条件是指在t=0时，y的值。

在MATLAB中，我们可以使用一个向量来表示初值条件。

例如，如果我们要求解以下常微分方程组：dy1/dt = -y1 + 2*y2dy2/dt = -2*y1 + 3*y2初值条件为：y(0) = [1; 0]我们可以定义初值条件为：y0 = [1; 0];现在，我们可以使用ode45函数来求解常微分方程组。

ode45函数的语法如下：[t,y] = ode45(f,tspan,y0)其中，f是一个函数句柄，tspan是一个包含起始时间和结束时间的向量，y0是一个包含初值条件的向量。

ode45函数将返回一个包含时间和解向量的矩阵。

例如，如果我们要求解以下常微分方程组：dy1/dt = -y1 + 2*y2dy2/dt = -2*y1 + 3*y2初值条件为：y(0) = [1; 0]时间范围为0到10秒，我们可以使用以下代码来求解：f = @(t,y) [-y(1) + 2*y(2); -2*y(1) + 3*y(2)];tspan = [0 10];y0 = [1; 0];[t,y] = ode45(f,tspan,y0);现在，我们可以绘制解向量随时间变化的图像。

matlab四阶龙格库塔法解方程组摘要：一、引言二、龙格库塔法介绍1.龙格库塔法的基本原理2.龙格库塔法的发展历程三、MATLAB 实现四阶龙格库塔法1.MATLAB 中龙格库塔法的函数2.四阶龙格库塔法的MATLAB 实现四、龙格库塔法解方程组的应用1.线性方程组的求解2.非线性方程组的求解五、结论正文：一、引言在数学领域，求解方程组是一项基本任务。

龙格库塔法作为高效数值求解线性方程组的方法，被广泛应用于各个领域。

MATLAB 作为一款强大的数学软件，可以方便地实现龙格库塔法求解方程组。

本文将介绍MATLAB 中四阶龙格库塔法解方程组的原理、实现与应用。

二、龙格库塔法介绍1.龙格库塔法的基本原理龙格库塔法（Runge-Kutta method）是一种求解常微分方程初值问题的数值方法。

它通过求解一组线性方程来逼近微分方程的解，具有较高的数值稳定性和精度。

龙格库塔法可以分为四阶、五阶等多种形式，其中四阶龙格库塔法是较为常用的一种。

2.龙格库塔法的发展历程龙格库塔法由德国数学家卡尔·龙格（Carl Runge）和英国数学家詹姆斯·库塔（James Kutta）分别在1900 年和1901 年独立发现。

自那时以来，龙格库塔法在数学、物理、工程等领域得到了广泛应用，并发展出了多种改进和扩展。

三、MATLAB 实现四阶龙格库塔法1.MATLAB 中龙格库塔法的函数在MATLAB 中，可以使用内置函数ode45、ode23、ode113 等实现龙格库塔法求解常微分方程。

这些函数分别对应四阶、五阶和三阶龙格库塔法。

2.四阶龙格库塔法的MATLAB 实现以下是一个使用MATLAB 实现四阶龙格库塔法求解方程组的示例：```matlabfunction [x, status] = solve_system_with_ode45(A, B, x0)% 定义方程组func = @(t, x) A * x + B;% 初始条件x0 = [1; 2];% 时间区间tspan = [0, 10];% 求解[x, status] = ode45(func, tspan, x0);end```四、龙格库塔法解方程组的应用1.线性方程组的求解线性方程组在数学、物理、工程等领域具有广泛应用。

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

它是由卡尔·迭龙格（Carl Runge）和马丁·威尔黑尔姆·库塔（Wilhelm Kutta）在20世纪初提出的，该方法以两位数值分析家的名字来命名。

二、简单描述迭龙格-库塔法是通过数值逼近的方式，来计算常微分方程的近似解。

它是一种显式求解方法，适用于解非线性常微分方程和具有较大阶数的常微分方程。

三、数学原理迭龙格-库塔法主要是通过将微分方程转化为差分方程，利用数值解的方式来逼近微分方程的解。

它是一种显式方法，通过不断迭代得到下一个时间步的近似解。

四、matlab中的应用在matlab中，可以使用ode45函数来调用迭龙格-库塔法求解常微分方程。

ode45函数是matlab中集成的一个函数，通过调用ode45函数，可以直接求解常微分方程的数值解。

五、实例演示下面通过一个简单的例子来演示如何使用matlab中的ode45函数来求解常微分方程。

我们考虑一个简单的一阶常微分方程：dy/dt = -y初始条件为y(0) = 1。

在matlab中，可以通过以下代码来求解该微分方程：```定义微分方程的函数function dydt = myode(t, y)dydt = -y;调用ode45函数求解[t, y] = ode45(myode, [0, 5], 1);plot(t, y);```运行以上代码，即可得到微分方程的数值解，并通过绘图来展示解的变化。

六、总结迭龙格-库塔法是一种常用的数值解常微分方程的方法，它在matlab中有较为方便的调用方式。

通过ode45函数，可以快速求解常微分方程的数值解，并通过绘图来展示结果。

希望本篇文章对读者有所帮助，谢谢阅读。

七、应用场景和优势在实际应用中，迭龙格-库塔法广泛应用于各种科学和工程领域，如物理学、化学、生物学、经济学等。

龙格库塔法（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);```。

-0.5
-0.55
-0.6
-0.65
-0.7
-0.75
-0.8
-0.85
-0.9
-0.95
-1
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
time t0=0,tt=1
图3 给定新的初始数据，由函数xprim2定义的ODE解的图形
（d） 求解下面方程组并不难：
x x x x ì ' = - 0.1
在下面的初值问题中，有两个未知函数：x1(t)和x2(t),并用以下式子表达其微... 页码，1/11
Matlab关于微分方程的解法
MATLAB使用龙格-库塔-芬尔格（Runge-Kutta-Fehlberg）方法来解ODE问题。在有限点内计算求解。而 这些点的间距有解的本身来决定。当解比较平滑时，区间内使用的点数少一些，在解变化很快时，区间内应使 用较多的点。 为了得到更多的有关何时使用哪种解法和算法的信息，推荐使用helpdesk。所有求解方程通用的语法或句法在 命令集中头两行给出。时间间隔将以向量t=[t0,tt]给出。 命令ode23可以求解(2,3)阶的常微分方程组，函数ode45使用(4,5)阶的龙格-库塔-芬尔格方法。注意，在这种情 况下x’是x的微分不是x的转置。 在命令集中solver将被诸如ode45函数所取代。
0.6
0.55
0.5
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
time t0=0,tt=1
图1 由函数xprim1定义的ODE解的图形
（b） 解下面的ODE过程是等价的：
ïíìx' = x2
ïîx(0) = 1

Matlab求解教程：求解延迟微分方程使用MA TLAB求解延时微分方程的两种方法：DDE23和SimuLink有些不同点需要注意，否则结果会出现错误使用MA TLAB来求解延迟微分方程是在生物数学和化学计算求解中经常遇到的事，在其它领域也比较常见。

我所知道的，在MA TLAB中求解微分方程有三种方法：1.使用ode45（龙格-库塔法的一个变种）求解，这时用一个数组，记录y的延迟项，但是c的值要考虑步长，再代入方程就能实现延时效应；2.使用dde23求解常数延时方程、使用ddesd可以用来求解延迟与时间t有关的延迟微分方程；3.使用SimuLink建模求解，SimuLink是求解广义微分代数系统的通用工具，功能很强大，但是看惯了编程指令的人可能不大习惯，调试似乎也不太方便。

既然本文专门讨论求解延迟微分方程，就先介绍一下专用工具dde23，dde系列求解函数是由Southern Methodist University 的L.F. Shampine 和S. Thompson根据他们早期使用Fortran编写的Fortran 90 DDE Solver 移植到MATLAB上的，从MA TLAB6.5开始加入MA TLAB的官方发行版，根据薛定宇教授在其几本关于MA TLAB的著作中提到的，该函数返回的sol中结构体sol.x 和sol.y均为按行排列，与ode45等不同，不太规范（没办法，因为这个函数本来就不是Mathworks的官方作品），不过这一点已经不大可能得到改进了，因为L.F. Shampine 和S. Thompson 已经决定停止维护这个文件。

如果您想进一步了解该函数，可以访问它的主页。

MA TLAB帮助中关于该函数的介绍不很清楚，如果需要进一步了解这个函数，需要下载作者为其写的手册。

下面以MA TLAB 中所附的一个例程来说明这个函数与Simulink建模求解的不同。

在MA TLAB prompt中键入edit ddex1就会找看到函数作者所写的一个入门例子：function ddex1%DDEX1 Example 1 for DDE23.% This is a simple example of Wille' and Baker that illustrates the% straightforward formulation, computation, and plotting of the solution% of a system of delay differential equations (DDEs).%% The differential equations%% y'_1(t) = y_1(t-1)% y'_2(t) = y_1(t-1)+y_2(t-0.2)% y'_3(t) = y_2(t)%% are solved on [0, 5] with history y_1(t) = 1, y_2(t) = 1, y_3(t) = 1 for% t <= 0.%% The lags are specified as a vector [1, 0.2], the delay differential% equations are coded in the subfunction DDEX1DE, and the history is% evaluated by the function DDEX1HIST. Because the history is constant it% could be supplied as a vector:% sol = dde23(@ddex1de,[1, 0.2],ones(3,1),[0, 5]);%% See also DDE23, FUNCTION_HANDLE.% Jacek Kierzenka, Lawrence F. Shampine and Skip Thompson% Copyright 1984-2004 The MathW orks, Inc.% $Revision: 1.2.4.2 $ $Date: 2005/06/21 19:24:16 $sol = dde23(@ddex1de,[1, 0.2],@ddex1hist,[0, 5]);figure;plot(sol.x,sol.y)title('An example of Wille'' and Baker.');xlabel('time t');ylabel('solution y');% --------------------------------------------------------------------------function s = ddex1hist(t)% Constant history function for DDEX1.s = ones(3,1);% --------------------------------------------------------------------------function dydt = ddex1de(t,y,Z)% Differential equations function for DDEX1.ylag1 = Z(:,1);ylag2 = Z(:,2);dydt = [ ylag1(1)ylag1(1) + ylag2(2)y(2) ];这里先不管函数使用的具体语法，求解模型为：显然有两个延时常数1、0.2。

1. matlab 新建.m文件，编写龙格-库塔法求解函数function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor((b-a)/h); %求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end2.另外再新建一个.,m文件，定义要求解的常微分方程函数function dx=fun1(t,x)dx =zeros(2,1);%初始化列向量dx(1) =0.08574*x(2)-1.8874*x(1)-20.17;dx(2) =1.8874*x(1)-0.08574*x(2)+115.16;3,再新建一个.m文件，利用龙格-库塔方法求解常微分方程[T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20); %求解步骤2定义的fun1常微分方程，@fun1是调用其函数指针，从0到20，间隔为1subplot(122)plot(T1,F1)%自编的龙格库塔函数效果title('自编的龙格库塔函数')grid运行步骤3文件即可得到结果，F1为估测值或者可以调用matlab自带函数ode45求解方法如下：[T,F] = ode45(@fun1,0:1:20,[17.64 37800 ]); subplot(121)plot(T,F)%Matlab自带的ode45函数效果title('ode45函数效果')grid。

matlab中三元二阶微分方程ode45编程在Matlab中，可以使用ode45函数来求解三元二阶微分方程。

ode45是一个常用的数值解微分方程的函数，它使用了一种基于龙格-库塔方法的算法来求解微分方程。

下面是使用ode45函数求解三元二阶微分方程的步骤：1. 定义一个匿名函数，该函数输入参数为t和y，其中t表示时间，y 表示状态变量向量。

该匿名函数需要返回状态变量向量y的一阶导数dy/dt和二阶导数d2y/dt2。

例如，假设要求解以下三元二阶微分方程：d2x/dt2 + 2*dx/dt + 3*x = sin(t)则可以定义如下匿名函数：f = @(t,y) [y(2); -2*y(2)-3*y(1)+sin(t)];其中，y(1)表示状态变量x，y(2)表示状态变量dx/dt。

2. 定义初始条件。

即定义状态变量向量在初始时刻的值。

例如：y0 = [0;0];表示在初始时刻，状态变量x和dx/dt均为0。

3. 定义时间区间。

即定义需要求解微分方程的时间区间。

例如：tspan = [0,10];表示需要求解从时间0到时间10的微分方程。

4. 调用ode45函数进行求解。

例如：[t,y] = ode45(f,tspan,y0);其中，t表示求解得到的时间向量，y表示求解得到的状态变量向量矩阵。

矩阵y的每一行对应一个时间点的状态变量值。

5. 绘制结果。

可以使用plot函数绘制状态变量随时间的变化曲线。

例如：plot(t,y(:,1),'-o',t,y(:,2),'-*');表示绘制状态变量x和dx/dt随时间的变化曲线。

以上就是使用ode45函数求解三元二阶微分方程的基本步骤。

需要注意的是，当定义匿名函数时，需要按照正确的语法规则编写微分方程，并确保输入参数和输出结果符合要求。

此外，在定义初始条件和时间区间时也需要注意输入参数的格式和范围。

总之，使用ode45函数求解三元二阶微分方程是一种简单而有效的数值解法，可以帮助工程师和科学家快速得到微分方程的数值解，并用于实际应用中。

(完整word 版)龙格库塔法求微分方程matlab龙格—库塔方法求解微分方程初值问题（数学1201+41262022+陈晓云）初值问题：y x x -+=2dxdy ，10≤≤x 1)0(y = 四阶龙格-库塔公式：()y x K n n ,f 1=⎪⎪⎪⎪⎭⎫ ⎝⎛+=+K h y x K n h n 122f ,2 ⎪⎭⎫ ⎝⎛++=K y x fK h n h n 232,2 ()K h y h x f K n n 34,++=()K K K K y y h n 43211n 226++++=+ 程序：1）建立四阶龙格-库塔函数function [ x,y ] = nark4( dyfun,xspan,y0,h )% dyfun 为一阶微分方程的函数；y0为初始条件；xspan 表示x 的区间；h 为区间的步长； x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);(完整word版)龙格库塔法求微分方程matlab k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+k2*2+2*k3+k4)/6;endx=x;y=y;2）执行程序（m文件）dyfun=inline('x^2+x-y');[x,y1]=nark4(dyfun,[0,1],1,0.1);x=0:0.1:1;Format longy2=x.^2-x+1R4=y2-y1[x',y1',y2',R4']y2=dsolve('Dy=x^2+x-y','y(0)=1','x')plot(x,y1,'b*-')hold ony3=inline('x^2-x+1')fplot(y3,[0,1],'ro-')legend('R-K4','解析解')3）执行结果ans =X RK4近似值解析值0 1.000000000000000 1.000000000000000 0.100000000000000 0.910000208333333 0.910000000000000 0.200000000000000 0.840000396841146 0.840000000000000 0.300000000000000 0.790000567410084 0.790000000000000 0.400000000000000 0.760000721747255 0.760000000000000 0.500000000000000 0.750000861397315 0.750000000000000 0.600000000000000 0.760000987757926 0.760000000000000 0.700000000000000 0.790001102093746 0.790000000000000 0.800000000000000 0.840001205549083 0.8400000000000000.900000000000000 0.910001299159352 0.9100000000000001.000000000000000 1.000001383861433 1.000000000000000误差-0.000000208333333-0.000000396841146-0.000000567410084-0.000000721747255-0.000000861397315-0.000000987757926-0.000001102093746-0.000001205549083-0.000001299159352-0.000001383861433 y2 =x^2 - x + 1结果分析：初值问题的解析解为Y=x^2 - x + 1由图看出龙格库塔方法误差很小，具有很高的精度。

matlab中常用来解常微分方程的原理Matlab是一种常用的科学计算软件，它在解常微分方程方面具有很强的功能。

在实际应用中，常微分方程是描述自然现象和工程问题的重要数学模型。

通过使用Matlab，我们可以方便地求解各种类型的常微分方程，从而得到问题的解析解或数值解。

本文将介绍Matlab中常用来解常微分方程的原理和方法。

常微分方程是描述一个未知函数及其导数之间关系的方程。

它可以分为初值问题和边值问题两种类型。

初值问题是给定方程的初始条件，要求求解出函数的解析解或数值解；而边值问题是给定方程在两个或多个点上的边界条件，要求求解出满足这些条件的函数解。

Matlab提供了多种方法来解这些问题，如常微分方程求解器、符号计算工具箱和数值计算工具箱等。

我们来介绍常微分方程求解器。

Matlab中的常微分方程求解器可以直接求解一阶或高阶的常微分方程。

常见的求解器有ode45、ode23、ode15s等。

其中，ode45是一种常用的求解器，它采用的是龙格-库塔方法，可以求解刚性和非刚性的常微分方程。

使用这些求解器，我们只需要定义好方程的形式和初始条件，就可以得到方程的数值解。

我们介绍符号计算工具箱。

符号计算工具箱可以对常微分方程进行符号计算，得到方程的解析解。

通过符号计算，我们可以得到方程的精确解，而不仅仅是数值解。

在Matlab中，符号计算工具箱提供了dsolve函数来求解常微分方程。

我们只需要输入方程的形式，就可以得到方程的解析解。

符号计算工具箱对于一些简单的常微分方程求解非常方便，但对于复杂的方程，可能需要更复杂的方法来求解。

我们介绍数值计算工具箱。

数值计算工具箱提供了各种数值方法来求解常微分方程。

常见的数值方法有欧拉方法、龙格-库塔方法和有限差分法等。

这些方法可以将常微分方程转化为一系列迭代计算，从而得到方程的数值解。

在Matlab中，数值计算工具箱提供了多个函数来实现这些数值方法，如euler、ode23s和fdcoeffF等。

ylabel('角度');
A=[x,y];
%y(find(y==max(y)))
%Num=(find(y==max(y)))
[y,T]=max(y);
fprintf('角度峰值等于%d',y) %角度的峰值也就是 π
fprintf('\n')
fprintf('周期等于%d',T*h)
legend('欧拉');
稳定性很好，RK4 法是四阶方法，每步的误差是 h5 阶，而总积累误差为 h4 阶。所以比欧 拉稳定。
运行第三个程序：在一幅图中显示欧拉法和 RK4 法，随着截断误差的积累，欧拉法产生了 较大的误差 h=0.01
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术关，系电，力根通保据过护生管高产线中工敷资艺设料高技试中术卷资，配料不置试仅技卷可术要以是求解指，决机对吊组电顶在气层进设配行备置继进不电行规保空范护载高与中带资负料荷试下卷高问总中题体资，配料而置试且时卷可，调保需控障要试各在验类最；管大对路限设习度备题内进到来行位确调。保整在机使管组其路高在敷中正设资常过料工程试况中卷下，安与要全过加，度强并工看且作护尽下关可都于能可管地以路缩正高小常中故工资障作料高；试中对卷资于连料继接试电管卷保口破护处坏进理范行高围整中，核资或对料者定试对值卷某，弯些审扁异核度常与固高校定中对盒资图位料纸置试，.卷保编工护写况层复进防杂行腐设自跨备动接与处地装理线置，弯高尤曲中其半资要径料避标试免高卷错等调误，试高要方中求案资技，料术编试交写5、卷底重电保。要气护管设设装线备备置敷4高、调动设中电试作技资气高，术料课中并3中试、件资且包卷管中料拒含试路调试绝线验敷试卷动槽方设技作、案技术，管以术来架及避等系免多统不项启必方动要式方高，案中为；资解对料决整试高套卷中启突语动然文过停电程机气中。课高因件中此中资，管料电壁试力薄卷高、电中接气资口设料不备试严进卷等行保问调护题试装，工置合作调理并试利且技用进术管行，线过要敷关求设运电技行力术高保。中护线资装缆料置敷试做设卷到原技准则术确：指灵在导活分。。线对对盒于于处调差，试动当过保不程护同中装电高置压中高回资中路料资交试料叉卷试时技卷，术调应问试采题技用，术金作是属为指隔调发板试电进人机行员一隔，变开需压处要器理在组；事在同前发一掌生线握内槽图部内 纸故，资障强料时电、，回设需路备要须制进同造行时厂外切家部断出电习具源题高高电中中源资资，料料线试试缆卷卷敷试切设验除完报从毕告而，与采要相用进关高行技中检术资查资料和料试检，卷测并主处且要理了保。解护现装场置设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。

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

matlab编的4阶龙格库塔法解微分方程的程序2010-03—10 20：16function varargout=saxplaxliu(varargin)clc，clearx0=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:nfprintf('%2d %4.4f ％4.4f\n’,i，x(i），y（i))；endfunction 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-1K1=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）;endy=y1；结果：i x（i） y（i)1 0.0000 1。

00002 0。

1000 0.99013 0。

2000 0.96154 0.3000 0。

91745 0.4000 0。

86216 0.5000 0。

80007 0。

6000 0。

73538 0。

7000 0.67119 0。

8000 0。

609810 0。

9000 0。

552511 1.0000 0.500012 1.1000 0.452513 1.2000 0.4098yi+1=yi+h＊（ K1+ K2)/2K1=f（xi,yi)K2=f(xi+h，yi+h*K1)依次类推，如果在区间［xi，xi+1］内多预估几个点上的斜率值K1、K2、……Km,并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。