微分方程的MATLAB实现与欧拉算法

数学实验—实验报告一、实验项目： 二、实验目的和要求1、本章将对人口变化、动物种群变迁、网络系统的可靠性分析，介绍微分方程（组）的模型建立、数值解和图形解等方法，并用MATLAB 几何直观地展示各种求解方法的求解结果。

2、利用欧拉公式求解方程三、实验题目问题一：求微分方程的解析解，并画出它们的图形， y ’=y +2x , y (0)=1, 0<x <1;问题二：用向前欧拉公式和改进的欧拉公式求方程y ’=y -2x /y , y (0)=1的数值解（0≤x ≤1 , h =0.1） 要求编写程序。

问题三：Rossler 微分方程组当固定参数b=2, c=4时，试讨论随参数a 由小到大变化（如a ∈(0,0.65]）而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？问题四：水的流出时间一横截面积为常数A ，高为H 的水池内盛满水，由池底一横截面积为B 的小孔放水。

设水从小孔流出的速度为v=(2gh)0.5，求在任意时刻的水面高度和将水放空所需的时间。

时间t →高度h 。

问题五：考虑相互竞争模型两种相似的群体之间为了争夺有限的同一种事物来源和生存空间而进行生死存亡竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从Logistics 规律。

三、实验过程问题一：用matlab 编写代码： x=[0,1]y=dsolve('Dy=y+2*x')y=dsolve('Dy=y+2*x', 'y(0)=1', 'x')ezplot(x,y)输出：y =-2*x+exp(t)*C1 （通解）y =-2*x-2+3*exp(x) 画图：x=0:0.01:1;y =-2*x-2+3*exp(x);plot(x,y)'''()x y z y x ayz b z x c =--⎧⎪=+⎨⎪=+-⎩0.10.20.30.40.50.60.70.80.91问题二： 1、分析：解：（1）解析解法得到其精确解：（2） 向前欧拉法：1(2/)n n n n n y y h y x y +=+- (1)2/n n nh y h x y =+- 迭代公式为 n+1y 1.10.2/n n n y x y =-，其中0y =y(0)=1（3）改进欧拉法：n+1nn n n n+1n +1n n n n n n n n n n n nn2n nnnn y =y +(h /2)*[(y -2x /y )+(y -2x /y )]=y +(h /2)*[(y -2x/y )+(y +h -2(x+h )/(yh ))] =y+(h /2)*[2y2x /y2(x +h )/(y+h )]=(1+h )y /2x/y(x +h )/(y+h )h h h h +--+-- 迭代公式为 n+1y 1.10.1/0.1()/()0.005n n n n n y x y x h y h =--+++，其中0y =y(0)=12、Matlab 编码x1(1)=0;y1(1)=1;y2(1)=1;h=0.1; for k=1:10 x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)+2*h*x1(k)/y1(k);y2(k+1)=(1+h)*y2(k)+(h*h)/2-h*x1(k)/y2(k)-h*(x1(k)+h)/(y2(k)+h); end x=0:0.1:1; y=(2*x+1).^(1/2);x1=x1(1:11),y=y(1:11),y1=y1(1:11),y2=y2(1:11), plot(x,y,x1,y1,'k:',x1,y2,'r--')显示图像及结果：x1 = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000y = 1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321y1 = 1.0000 0.9000 0.8322 0.7971 0.7926 0.8143 0.8557 0.9103 0.9731 1.0402 1.1092y2 = 1.0000 1.0959 1.1847 1.2679 1.3468 1.4222 1.4948 1.5653 1.6340 1.7016 1.768300.10.20.30.40.50.60.70.80.910.811.21.41.61.82图中，蓝色曲线是精确解，红色曲线是向前欧拉法曲线，黑色曲线是改进后欧拉法曲线问题三1、matlab 编程function r=rossler(t,x) global a; global b; global c;r=[-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c)];global a; global b; global c; b=2; c=4; t0=[0,200]; for a=0:0.02:0.65[t,x]=ode45('rossler',t0,[0,0,0]); subplot(1,2,1);plot(t,x(:,1),'r',t,x(:,2),'g',t,x(:,3),'b');title('x(ºìÉ«),y(ÂÌÉ«),z(ÀºÉ«)Ëæt±ä»¯Çé¿ö');xlabel('t'); subplot(1,2,2);plot3(x(:,1),x(:,2),x(:,3))title('Ïàͼ');xlabel('x');ylabel('y');zlabel('z'); pause end当a=0时，图像如下50100150200-1-0.8-0.6-0.4-0.200.20.40.6x(红色),y(绿色),z(篮色)随t 变化情况t-0.5相图z当a=0时，(x,y,z)收敛于（0,0.5,0.5）当a=0.12时，图像如下50100150200-1.2-1-0.8-0.6-0.4-0.200.20.40.6x(红色),y(绿色),z(篮色)随t 变化情况t-0.5相图z当a=0.28时，(x,y ,z)仍然收敛，但是收敛速度大大降低。

为了更好地理解欧拉法求解一阶微分方程在Matlab中的应用，我们首先来了解一些背景知识。

一阶微分方程是指只含有一阶导数的方程，通常表示为dy/dx=f(x,y)，其中f(x,y)是关于x和y的函数。

欧拉法是一种常见的数值解法，用于求解微分方程的近似数值解。

它是一种基本的显式数值积分方法，通过将微分方程转化为差分方程来进行逼近。

在Matlab中，我们可以利用欧拉法求解一阶微分方程。

我们需要定义微分方程的函数表达式，然后选择合适的步长和初始条件，最后使用循环计算逼近解。

下面我们来具体讨论如何在Matlab中使用欧拉法来求解一阶微分方程。

我们假设要求解的微分方程为dy/dx=-2x+y，初始条件为y(0)=1。

我们可以通过以下步骤来实现：1. 我们需要在Matlab中定义微分方程的函数表达式。

在Matlab中，我们可以使用function关键字来定义函数。

在这个例子中，我们可以定义一个名为diff_eqn的函数，表示微分方程的右侧表达式。

在Matlab中，这个函数可以定义为：```matlabfunction dydx = diff_eqn(x, y)dydx = -2*x + y;end```2. 我们需要选择合适的步长和初始条件。

在欧拉法中，步长的选择对于数值解的精度非常重要。

通常情况下，可以先尝试较小的步长，然后根据需要进行调整。

在这个例子中，我们可以选择步长h=0.1，并设置初始条件x0=0，y0=1。

3. 接下来，我们可以使用循环来逼近微分方程的数值解。

在每一步，根据欧拉法的迭代公式y(i+1) = y(i) + h * f(x(i), y(i))，我们可以按照下面的Matlab代码计算逼近解：```matlabh = 0.1; % 步长x = 0:h:2; % 定义计算区间y = zeros(1, length(x)); % 初始化y的值y(1) = 1; % 设置初始条件for i = 1:(length(x)-1) % 欧拉法迭代y(i+1) = y(i) + h * diff_eqn(x(i), y(i));end```通过上述步骤，在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求微分方程数值解利用matlab求微分方程数值解是一种常用的数学计算方法。

在实际工程和科学研究中，许多问题都可以用微分方程来描述，但是解析解往往难以求得，因此需要用数值方法求解微分方程。

求解微分方程的数值方法有很多种，其中比较常用的是欧拉法和龙格-库塔法。

欧拉法是一种基本的数值方法，它采用离散化的方法将微分方程转化为差分方程，然后通过迭代来求出数值解。

欧拉法的具体步骤是：首先将自变量和因变量离散化，然后利用微分方程的定义式将微分方程转化为差分方程，最后通过迭代求出数值解。

欧拉法的优点是简单易懂，但是精度较低，容易产生误差。

龙格-库塔法是一种高阶数值方法，它将微分方程转化为一系列的差分方程，并采用递推的方法求解数值解。

龙格-库塔法的优点是精度高，收敛速度快，适用于求解复杂的微分方程。

但是龙格-库塔法的计算量较大，需要进行多次计算，计算时间较长。

在使用matlab求解微分方程时，可以直接调用matlab中的ode 函数来求解微分方程。

ode函数是matlab中内置的求解微分方程的函数，它支持多种数值方法，包括欧拉法和龙格-库塔法等。

使用ode函数可以简化求解微分方程的过程，提高计算效率。

在使用ode函数求解微分方程时，需要先定义微分方程的函数表达式，然后将函数表达式作为参数传入ode函数中。

ode函数会自动选择合适的数值方法来求解微分方程，并返回数值解。

通过调整ode函数的参数，可以进一步提高求解微分方程的精度和计算效率。

除了ode函数外，matlab中还有很多其他的数值计算函数，如dsolve函数、pdepe函数等，它们可以用来求解不同类型的微分方程。

在实际应用中，需要根据具体问题选择合适的数值方法和函数来求解微分方程。

利用matlab求解微分方程数值解是一种常用的数学计算方法，可以通过调用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);```。

euler方法求解常微分方程matlab以Euler方法求解常微分方程(Matlab)在数学和工程领域中，常微分方程是一类描述物体运动、电路行为、传热传质等现象的数学模型。

求解常微分方程有许多方法，其中一种简单而常用的方法是Euler方法。

本文将介绍如何使用Matlab编程语言来实现Euler方法求解常微分方程，并给出一个具体的例子进行说明。

Euler方法基本思想是将微分方程转化为差分方程，通过一系列的逼近来得到方程的近似解。

对于一阶常微分方程dy/dx=f(x,y)，假设在区间[a,b]上等距离取n个点，则每个点的步长为h=(b-a)/n。

对于每个点(xi,yi)，根据微分方程可以得到差分方程yi+1 = yi + h * f(xi,yi)，其中f(xi,yi)代表在点(xi,yi)处的斜率。

下面我们以一个简单的例子来说明Euler方法的应用。

考虑一阶线性微分方程dy/dx = -2xy，初始条件为y(0)=1，求解在区间[0,1]上的近似解。

我们需要将微分方程转化为差分方程。

根据Euler方法的思想，我们有yi+1 = yi + h * (-2xi * yi)。

然后，我们需要选择合适的步长h和区间[0,1]上的点数n。

在Matlab中，我们可以定义一个函数来表示微分方程。

代码如下：function dydx = myODE(x,y)dydx = -2*x*y;end然后，我们可以编写一个求解Euler方法的函数。

代码如下：function [x,y] = myEuler(a,b,n,y0)h = (b-a)/n;x = a:h:b;y = zeros(1,n+1);y(1) = y0;for i = 1:ny(i+1) = y(i) + h * myODE(x(i),y(i));endend我们可以调用myEuler函数来求解微分方程。

代码如下：a = 0;b = 1;n = 100;y0 = 1;[x,y] = myEuler(a,b,n,y0);这样，我们就得到了在区间[0,1]上的近似解。

欧拉法是数值分析中常用的一种方法，用于求解常微分方程的数值解。

在MATLAB中，可以通过编写相应的代码来实现欧拉法求解微分方程。

下面我们将通过具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。

我们要了解欧拉法的基本原理。

欧拉法是一种通过迭代逼近微分方程解的方法，它基于微分方程的定义，通过离散化的方法逼近微分方程的解。

其基本思想是利用微分方程的导数定义，将微分方程以差分形式进行逼近。

具体而言，欧拉法通过将微分方程转化为差分方程的形式，然后通过迭代逼近得到微分方程的数值解。

接下来，我们通过一个具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。

假设我们要求解以下的一阶常微分方程：(1) dy/dx = x + y(2) y(0) = 1现在我们来编写MATLAB代码来实现欧拉法求解这个微分方程。

我们需要确定微分方程的迭代步长和迭代范围。

假设我们将x的范围取为0到10，步长为0.1。

接下来，我们可以编写MATLAB代码如下：```matlab欧拉法求解微分方程 dy/dx = x + y定义迭代步长和范围h = 0.1;x = 0:h:10;初始化y值y = zeros(1,length(x));y(1) = 1;使用欧拉法迭代求解for i = 1:(length(x)-1)y(i+1) = y(i) + h * (x(i) + y(i));end绘制图像plot(x,y,'-o');xlabel('x');ylabel('y');title('欧拉法求解微分方程 dy/dx = x + y');```在这段MATLAB代码中，我们首先定义了迭代的步长和范围，并初始化了微分方程的初始值y(0) = 1。

然后通过for循环使用欧拉法进行迭代求解微分方程，最后绘制出了微分方程的数值解的图像。

通过以上的实例讲解，我们可以看到，在MATLAB中使用欧拉法求解微分方程是非常简单而直观的。

Matlab中欧拉法求解常微分方程初值问题一、概念介绍在数学和工程领域，常微分方程初值问题是一个广泛应用的数学概念。

它描述了一个未知函数在给定初始条件下的行为。

而欧拉法则是一种常用的数值方法，用来解决常微分方程初值问题。

在Matlab中，我们可以利用欧拉法来求解常微分方程问题，从而得到函数在给定初始条件下的近似解。

二、欧拉法的基本原理欧拉法的基本思想是通过离散化微分方程，将其转化为递推的差分方程。

考虑一个一阶常微分方程初值问题：\[ \frac{dy}{dx} = f(x, y), \quad y(x_0) = y_0 \]在欧拉法中，我们采用递推的方式，根据已知的初始条件和微分方程的性质，通过迭代来得到逼近解的数值结果。

具体地，我们首先将自变量$x$的范围进行等间距分割，得到$x_0, x_1, x_2, ..., x_n$，并将步长记为$h$。

根据微分方程的性质，我们可以根据已知的初始条件$y(x_0) = y_0$，通过迭代计算得到近似解$y(x_1), y(x_2), ..., y(x_n)$。

三、Matlab中的欧拉法求解在Matlab中，我们可以利用欧拉法来求解常微分方程初值问题。

以求解一阶常微分方程为例，假设我们需要求解以下的常微分方程初值问题：\[ \frac{dy}{dx} = -2xy, \quad y(0) = 1 \]我们可以利用欧拉法的思想，将自变量$x$的范围进行离散化，然后根据欧拉法的递推公式，利用迭代的方式得到近似解的数值结果。

具体地，在Matlab中，我们可以编写如下代码来实现欧拉法的求解过程：```matlabfunction y = euler_method(f, x0, y0, h, n)% 初始化存储结果的数组x = zeros(1, n+1);y = zeros(1, n+1);% 将初始条件存入数组x(1) = x0;y(1) = y0;% 利用欧拉法进行迭代for i = 1:nx(i+1) = x(i) + h;y(i+1) = y(i) + h * f(x(i), y(i));end% 返回近似解的数值结果plot(x, y); % 绘制解的图像end```在上述代码中，我们定义了一个名为`euler_method`的函数，其中包含了欧拉法的计算过程。

一、概述Matlab作为一种常用的科学计算软件，在微分方程的数值解法领域具有广泛的应用。

微分方程是描述自然现象中变化规律的数学工具，而数值解法则是指使用计算机进行近似求解微分方程的方法。

在Matlab 中，有多种常用的数值解法可以用来求解微分方程，例如欧拉法、改进的欧拉法、四阶龙格-库塔法等。

本文将对这些数值解法进行介绍和比较，以帮助读者更好地理解和应用微分方程求解数值方法。

二、欧拉法欧拉法是微分方程的最简单的数值解法之一，它通过离散化微分方程进行近似求解。

具体而言，对于一阶常微分方程dy/dx=f(x,y)，可以利用欧拉法进行数值解。

欧拉法的基本思想是将自变量x的增量Δx分成n个小区间，然后根据微分方程的数值近似公式y(x+Δx)=y(x)+f(x,y)Δx对每个小区间进行迭代计算。

欧拉法的优点是简单易实现，但由于它是一阶的数值方法，因此对于某些微分方程求解效果可能不够准确。

三、改进的欧拉法改进的欧拉法是对欧拉法的一种改进，它通过在每个小区间内使用平均斜率来提高求解的精度。

具体而言，对于微分方程dy/dx=f(x,y)，改进的欧拉法可以通过以下迭代公式进行数值求解：y(x+Δx)=y(x)+Δx/2[f(x,y)+f(x+Δx,y+Δx*f(x,y))]改进的欧拉法相比于欧拉法具有更高的数值精度，但计算量也相对增加。

四、四阶龙格-库塔法四阶龙格-库塔法是一种常用的数值微分方程求解方法，它通过四次迭代计算来获得微分方程的数值解。

具体而言，对于微分方程dy/dx=f(x,y)，四阶龙格-库塔法可以用以下公式进行数值求解：k1=f(x,y)k2=f(x+Δx/2,y+Δx/2*k1)k3=f(x+Δx/2,y+Δx/2*k2)k4=f(x+Δx,y+Δx*k3)y(x+Δx)=y(x)+Δx/6*(k1+2*k2+2*k3+k4)四阶龙格-库塔法相比于欧拉法和改进的欧拉法具有更高的数值精度和稳定性，但计算量也相对较大。

解微分方程的欧拉法，龙格-库塔法及其MATLAB简单实例欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前进EULER法、后退EULER法、改进的EULER法。

缺点：欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。

因此欧拉格式一般不用于实际计算。

改进欧拉格式：为提高精度，需要在欧拉格式的基础上进行改进。

采用区间两端的斜率的平均值作为直线方程的斜率。

改进欧拉法的精度为二阶。

算法为：微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。

对于常微分方程：?x∈[a,b]y(a) = y0可以将区间[a,b]分成n段，那么方程在第xi点有y'(xi) = f(xi,y(xi))，再用向前差商近似代替导数则为：在这里，h是步长，即相邻两个结点间的距离。

因此可以根据xi点和yi点的数值计算出yi+1来：i=0,1,2,L这就是向前欧拉格式。

改进的欧拉公式：将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均，即上式便是梯形的欧拉公式。

可见，上式是隐式格式，需要迭代求解。

为了便于求解，使用改进的欧拉公式：数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。

龙格-库塔方法的基本思想：在区间[xn,xn+1]内多取几个点，将他们的斜率加权平均，作为导数的近似。

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

令初值问题表述如下。

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

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

解微分方程的欧拉法，龙格-库塔法及其MATLAB简单实例欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前进EULER法、后退EULER法、改进的EULER法。

缺点：欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。

因此欧拉格式一般不用于实际计算。

改进欧拉格式：为提高精度，需要在欧拉格式的基础上进行改进。

采用区间两端的斜率的平均值作为直线方程的斜率。

改进欧拉法的精度为二阶。

算法为：微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。

对于常微分方程：?x∈[a,b]y(a) = y0可以将区间[a,b]分成n段，那么方程在第xi点有y'(xi) = f(xi,y(xi))，再用向前差商近似代替导数则为：在这里，h是步长，即相邻两个结点间的距离。

因此可以根据xi点和yi点的数值计算出yi+1来：i=0,1,2,L这就是向前欧拉格式。

改进的欧拉公式：将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均，即上式便是梯形的欧拉公式。

可见，上式是隐式格式，需要迭代求解。

为了便于求解，使用改进的欧拉公式：数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。

龙格-库塔方法的基本思想：在区间[xn,xn+1]内多取几个点，将他们的斜率加权平均，作为导数的近似。

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

令初值问题表述如下。

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

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

欧拉法(euler)求解常微分方程的matlab程序及案例欧拉法是一种常见的求解常微分方程的数值解法，在MATLAB中可以通过编写简单的程序实现。

本文将介绍欧拉法的MATLAB程序及应用案例。

首先，让我们考虑以下的常微分方程：dy/dx = f(x, y)其中y是关于x的函数，f是已知的函数。

我们可以通过欧拉法求解该方程。

欧拉法的基本思想是将区间[x0, xn]分成n等份，然后用以下式子计算y的值：y(i+1) = y(i) + h*f(x(i), y(i))其中h是步长，x(i)和y(i)分别表示当前的x和y值，y(i+1)表示下一个y值。

通过重复上述计算，欧拉法可以求出y在x=n处的值。

下面是欧拉法的MATLAB程序：% 默认参数x0 = 0; % 初始值xn = 1; % 终止值y0 = 1; % 初始y值h = 0.1; % 步长f = @(x, y) -y; % 函数n = (xn - x0) / h; % 时间步数x = x0; % 初始x值y = y0; % 初始y值for i = 1:ny = y + h * f(x, y);x = x + h;enddisp(['y在x = ', num2str(xn), '处的值为：',num2str(y)]);在上述程序中，我们定义了默认的初始值、终止值、初始y值和函数。

程序中的n表示时间步数，x和y分别表示当前的x和y值。

通过for循环，欧拉法可以重复计算y的值，并最终求出y在x=n处的值。

下面是一个用欧拉法求解dy/dx = -y的应用案例：% 默认参数x0 = 0; % 初始值xn = 5; % 终止值y0 = 1; % 初始y值h = 0.1; % 步长f = @(x, y) -y; % 函数n = (xn - x0) / h; % 时间步数x = x0; % 初始x值y = y0; % 初始y值% 初始化结果数组result = zeros(n + 1, 2);result(1,:) = [x0 y0];for i = 1:ny = y + h * f(x, y);x = x + h;% 保存结果result(i + 1,:) = [x y];end% 绘制图形plot(result(:,1), result(:,2), '-o');xlabel('x');ylabel('y');title('欧拉法求解dy/dx=-y');在上述案例中，我们使用默认的参数，求解dy/dx=-y的方程。

Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。

由角动量定理我们知道εJ M =化简得到0sin 22=+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。

实际上这是一个解二阶常微分方程的问题。

θ在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程W T =h mg ∆=2J 21ω化简得到四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。

所以比欧拉稳定。

运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。

通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。

三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');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('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%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)-1K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1 and L1K2=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;% K4L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);and L4RK_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);endplot(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)两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;%欧拉的初值Euler_z(1)=0;RK_y(1)=0;%龙格库塔初值RK_z(1)=0;for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 andL1K2=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);% K2 and L2K3=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);% K3 and L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); K4 and L4RK_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);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);%角度的峰值也就是πfprintf('角度峰值等于%d',y)fprintf('\n')%周期fprintf('周期等于%d',T*h)xlabel('时间');ylabel('角度');legend('欧拉','RK4');。

Euler法解常微分方程Euler法解常微分方程算法：Step 1 分别取积分上限、积分下限、步长Step 2计算判断是否成立，成立转到Step 3，否则继续进行Step 4Step 3 计算Step 4Euler法解常微分方程算程序：function euler2(fun,y0,A,h)%fun--y'%y0---初值%A----x取值范围%a----x左区间端点值%b----x右区间端点值%h----给定步长x=min(A);b=max(A);y=y0;while x<b-hb=y;y=y+h*feval(fun,x,b)x=x+h;end例：用Euler法计算下列初值问题（取步长h=0.2）输入：fun=inline('-y-x*y^2')euler2(fun,1,[0 0.6],0.2)得到：y =0.8000y =0.6144y =0.4613指导教师：年月日改进Euelr法解常微分方程改进Euler法解常微分方程算法：Step 1 分别取积分上限、积分下限、步长Step 2 取一个以h为步长，a，b分别为左右端点的矩阵Step 3 （1）做显性Euler预测（2）将带入Step 4计算判断是否成立，成立返回Step 3，否则继续进行Step 5 Step 5改进Euler法解常微分方程算程序：function gaijineuler2(fun,y0,A,h)%fun--y'%y0---初值%A----x取值范围%a----x左区间端点值%b----x右区间端点值%h----给定步长a=min(A);b=max(A);x=a:h:b;y(1)=y0;for i=1:length(x)-1w1=feval(fun,x(i),y(i));y(i+1)=y(i)+h*w1;w2=feval(fun,x(i+1),y(i+1));y(i+1)=y(i)+h*(w1+w2)/2;endx=x'y=y'例：用改进Euler法计算下列初值问题（取步长h=0.25）输入：fun=inline('-x*y^2')gaijineuler2(fun,2,[0 5],0.25)得到：x =0.25000.50000.75001.00001.25001.50001.75002.00002.25002.50002.75003.00003.25003.50003.75004.00004.25004.50004.75005.0000y =2.00001.87501.59391.28241.00960.79320.62820.50370.40970.33790.28240.23890.20430.17650.15380.13520.11960.10660.09550.08610.0779指导教师：年月日禁止复制北京石油化工学院Onlyunited。

数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。

由角动量定理我们知道εJ M =化简得到 0sin 22=+θθlg dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为θ，这样比较容易解。

实际上这是一个解二阶常微分方程的问题。

在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程 WT = h mg ?=2J 21ω 化简得到θθcos35499.722=dtd 重力加速度取9.806651使用欧拉法令dxdy z =，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前Euler 方法数值求解。

y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));y(0)=0z(0)=0精度随着h 的减小而更高，因为向前欧拉方法的整体截断误差与h 同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。

2.RK4-四阶龙格库塔方法使用四级四阶经典显式Rungkutta 公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。

所以比欧拉稳定。

运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。

通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。

三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');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('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%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)-1K1=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 L4RK_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);endplot(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) %周期两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;Euler_z(1)=0; %欧拉的初值RK_y(1)=0;RK_z(1)=0; %龙格库塔初值for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔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);% K2 and L2K3=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 );% K3 and L3K4=RK_z(i)+h*L3;L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4RK_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);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfpri ntf('\n') fprintf('周期等于%d',T*h) %周期xlabel('时间');ylabel('角度');legend('欧拉','RK4');。