MATLAB求解常微分方程数值解
Matlab中常用的数值计算方法
Matlab中常用的数值计算方法数值计算是现代科学和工程领域中的一个重要问题。
Matlab是一种用于数值计算和科学计算的高级编程语言和环境，具有强大的数值计算功能。
本文将介绍Matlab中常用的数值计算方法，包括数值积分、数值解微分方程、非线性方程求解和线性方程组求解等。
一、数值积分数值积分是通过数值方法来近似计算函数的定积分。
在Matlab中，常用的数值积分函数是'quad'和'quadl'。
'quad'函数可以用于计算定积分，而'quadl'函数可以用于计算无穷积分。
下面是一个使用'quad'函数计算定积分的例子。
假设我们想计算函数f(x) = x^2在区间[0, 1]上的定积分。
我们可以使用如下的Matlab代码：```f = @(x) x^2;integral = quad(f, 0, 1);disp(integral);```运行这段代码后，我们可以得到定积分的近似值，即1/3。
二、数值解微分方程微分方程是描述自然界各种变化规律的数学方程。
在科学研究和工程应用中，常常需要求解微分方程的数值解。
在Matlab中，可以使用'ode45'函数来求解常微分方程的数值解。
'ode45'函数是采用基于Runge-Kutta方法的一种数值解法。
下面是一个使用'ode45'函数求解常微分方程的例子。
假设我们想求解一阶常微分方程dy/dx = 2*x，初始条件为y(0) = 1。
我们可以使用如下的Matlab代码：```fun = @(x, y) 2*x;[x, y] = ode45(fun, [0, 1], 1);plot(x, y);```运行这段代码后，我们可以得到微分方程的数值解，并绘制其图像。
三、非线性方程求解非线性方程是指方程中包含非线性项的方程。
在很多实际问题中，我们需要求解非线性方程的根。
重要：MATLAB常微分方程(组)数值解法
Matlab常微分方程求解问题分类
边值问题：
初值问题：
• 定解附加条件在自变量 的一端
• 一般形式为： y' f (x, y)
y(a)
y0
• 初值问题的数值解法一 般采用步进法，如 Runge-Kutta法
➢ 在自变量两端均给定附加 条件
y' f (x, y)
➢ 一般形式：y(a)y1, y(b)y2
1.根据常微分方程要求的求解精度与速度要求
求解初值问题：
y
'
y
2x y
y ( 0 ) 1
(0x1)
比较ode45和ode23的求解精度和速度
ode45和ode23的比较-1
function xODE clear all clc
format long
y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y')
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances dCdt = [rA; rB; rC; rD; rE];
三个串联的CSTR等温反应器（例4-3）
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];
使用Matlab进行微分方程求解的方法
使用Matlab进行微分方程求解的方法引言微分方程是数学中非常重要的一部分，广泛应用于物理、经济、工程等领域。
对于大部分微分方程的解析解往往难以求得，而数值解法则成为了一种常用的解决手段。
Matlab作为一种强大的科学计算软件，也提供了丰富的工具和函数用于求解微分方程，本文将介绍一些常见的使用Matlab进行微分方程求解的方法。
一、数值求解方法1. 欧拉方法欧拉方法是最简单的一种数值求解微分方程的方法，它将微分方程的微分项用差分的方式进行近似。
具体的公式为：y(n+1) = y(n) + hf(x(n), y(n))其中，y(n)表示近似解在第n个点的值，h为步长，f(x, y)为微分方程的右端项。
在Matlab中使用欧拉方法进行求解可以使用ode113函数，通过设定不同的步长，可以得到不同精度的数值解。
2. 中点法中点法是较为精确的一种数值求解微分方程的方法，它的计算公式为：k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)y(n+1) = y(n) + k2中点法通过计算两个斜率的平均值来得到下一个点的值，相较于欧拉方法，中点法能提供更精确的数值解。
3. 4阶龙格库塔法龙格库塔法是一类高阶数值求解微分方程的方法，其中4阶龙格库塔法是最常用的一种。
它的计算公式为：k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)k3 = hf(x(n) + h/2, y(n) + k2/2)k4 = hf(x(n) + h, y(n) + k3)y(n+1) = y(n) + (k1 + 2k2 + 2k3 + k4)/64阶龙格库塔法通过计算多个斜率的加权平均值来得到下一个点的值，相较于欧拉方法和中点法，它的精度更高。
二、Matlab函数和工具除了可以使用以上的数值方法进行微分方程求解之外，Matlab还提供了一些相关的函数和工具，方便用户进行微分方程的建模和求解。
matlab 二阶常微分方程数值求解函数
matlab 二阶常微分方程数值求解函数【最新版】目录1.Matlab 二阶常微分方程数值求解函数概述2.二阶常微分方程的一般形式3.Matlab 中用于数值求解二阶常微分方程的函数4.数值求解的步骤5.结论正文Matlab 二阶常微分方程数值求解函数概述二阶常微分方程是指具有以下形式的微分方程：a * y"" + b * y" + c * y = f(x)。
其中，a、b、c 为常数，y 是函数，x 是自变量，f(x) 是已知函数。
求解这类微分方程对于许多实际问题具有重要意义，如物理、生物学、经济学等领域。
Matlab 作为一种广泛应用于科学计算的语言，提供了丰富的函数库用于数值求解二阶常微分方程。
二阶常微分方程的一般形式在 Matlab 中，二阶常微分方程的一般形式可以表示为：y"" + p(x) * y" + q(x) * y = r(x)其中，p(x)、q(x) 和 r(x) 是已知函数，y 是待求解的函数。
Matlab 中用于数值求解二阶常微分方程的函数Matlab 提供了多个函数用于数值求解二阶常微分方程，如 ode45、ode23、ode113 等。
这些函数的用法及参数如下：- ode45：该函数是四阶龙格库塔法（RK45）的实现，适用于大多数情况。
其用法为：解 = ode45(函数句柄，[a, b], [c, d], e, [f, g])其中，函数句柄是一个函数句柄，它接受自变量 x 和时间 t 作为参数，并返回因变量 y。
a 和 b 分别是函数 y 的第一个和第二个导数。
c 和 d 分别是函数 y 的第三个和第四个导数。
e 是初始条件。
f 和g 是边界条件。
- ode23：该函数是二阶龙格库塔法（RK23）的实现，适用于某些特殊情况。
其用法与 ode45 类似。
- ode113：该函数是十一阶龙格库塔法（RK113）的实现，适用于要求高精度解的情况。
matlab_常微分方程数值解法
dt 2
简朴问题可以求得解析解，多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn，用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题，因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明
型
ode45
非刚性 单步法；4，5 阶 R-K 措施；合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法；2，3 阶 R-K 措施；合计 使用于精度较低旳情形
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中对常微分方程进行数值求解的快速方法。
1. 基本概念在MATLAB中，可以使用ode45函数来对常微分方程进行数值求解。
ode45是一种常用的Runge-Kutta法，它可以自适应地选取步长，并且具有较高的数值精度。
使用ode45函数可以方便地对各种类型的常微分方程进行求解，包括一阶、高阶、常系数和变系数的微分方程。
2. 函数调用要使用ode45函数进行常微分方程的数值求解，需要按照以下格式进行函数调用：[t, y] = ode45(odefun, tspan, y0)其中，odefun表示用于描述微分方程的函数，tspan表示求解的时间跨度，y0表示初值条件，t和y分别表示求解得到的时间序列和对应的解向量。
3. 示例演示为了更好地理解如何使用ode45函数进行常微分方程的数值求解，下面我们以一个具体的例子来进行演示。
考虑如下的一阶常微分方程：dy/dt = -2*y其中，y(0) = 1。
我们可以编写一个描述微分方程的函数odefun：function dydt = odefun(t, y)dydt = -2*y;按照上述的函数调用格式，使用ode45函数进行求解：tspan = [0 10];y0 = 1;[t, y] = ode45(odefun, tspan, y0);绘制出解曲线：plot(t, y);4. 高级用法除了基本的函数调用方式外，MATLAB中还提供了更多高级的方法来对常微分方程进行数值求解。
可以通过设定选项参数来控制数值求解的精度和稳定性，并且还可以对刚性微分方程进行求解。
5. 性能优化在实际工程应用中，常常需要对大规模的常微分方程进行数值求解。
matlab 解常微分方程
matlab 解常微分方程Matlab是一种功能强大的数学软件，它提供了解常微分方程的工具和函数。
常微分方程是数学中的一种重要的方程类型，描述了各种物理、工程和生物现象的变化规律。
本文将介绍如何使用Matlab 解常微分方程，并通过具体的实例来说明其应用。
我们需要了解常微分方程的基本概念。
常微分方程是指一个函数的导数与自变量之间的关系方程。
常微分方程的解是该函数在给定初始条件下的解析解或数值解。
在Matlab中，我们可以使用ode45函数来求解常微分方程的数值解。
接下来，我们将以一个简单的一阶常微分方程为例来说明Matlab 的使用。
考虑以下的一阶常微分方程：dy/dx = x^2 - y我们将该方程转化为Matlab中的函数形式，并设定初始条件y(0) = 1。
代码如下：```matlabfunction dydx = myODE(x, y)dydx = x^2 - y;endxspan = [0 10];y0 = 1;[x, y] = ode45(@myODE, xspan, y0);plot(x, y)xlabel('x')ylabel('y')title('Solution of dy/dx = x^2 - y')```在上述代码中，我们首先定义了一个名为myODE的函数，该函数接受两个参数x和y，并返回dy/dx的值。
然后，我们使用ode45函数来求解该常微分方程的数值解。
最后，我们绘制了解的曲线图，并添加了相应的坐标轴标签和标题。
通过运行上述代码，我们可以得到常微分方程dy/dx = x^2 - y的数值解，并绘制出解的曲线图。
这个例子展示了Matlab解常微分方程的基本步骤和方法。
除了一阶常微分方程，Matlab还可以解决更高阶的常微分方程。
对于高阶常微分方程，我们可以将其转化为一组一阶常微分方程，并使用类似的方法来求解。
Matlab提供了一系列的函数和工具箱来处理不同类型的常微分方程，并提供了丰富的文档和示例来帮助用户理解和应用这些工具。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
利用MATLAB求解常微分方程数值解目录1. 内容简介 (1)2. Euler Method(欧拉法)求解 (1)2.1. 显式Euler法和隐式Euler法 (2)2.2. 梯形公式和改进Euler法 (3)2.3. Euler法实用性 (4)3. Runge-Kutta Method(龙格库塔法)求解 (5)3.1. Runge-Kutta基本原理 (5)3.2. MATLAB中使用Runge-Kutta法的函数 (7)4. 使用MATLAB求解常微分方程 (7)4.1. 使用ode45函数求解非刚性常微分方程 (8)4.2. 刚性常微分方程 (9)5. 总结 (9)参考文献 (11)附录 (12)1. 显式Euler法数值求解 (12)2. 改进Euler法数值求解 (12)3. 四阶四级Runge-Kutta法数值求解 (13)4.使用ode45求解 (14)1.内容简介把《高等工程数学》看了一遍，增加对数学内容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。
理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。
实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。
把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。
文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。
最后考察MATLAB中各个函数的适用范围，当遇到实际工程问题时能够正确地得到问题的数值解。
2.Euler Method(欧拉法)求解Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节内容分别介绍这3种方法的具体内容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。
本节考虑实际初值问题使用解析法，对方程两边同乘以得到下式两边同时求积分并采用分部积分得到解析解：本节后面将对此方程进行求解，并与精确解进行对比，分析Euler的可行性。
2.1.显式Euler法和隐式Euler法显式和隐式Euler法都属于一阶方法，显式Euler法的迭代公式简单，如下所示：对过上述公式对式进行迭代，其中步长，计算之间的数值，迭代求解的MATLAB程序见附录，能够得出精确解和数值解的图像，如图所示。
图2.1 显式Euler法精确解和数值解图像从图2.1中可以看出，显式Euler法在斜率很大的时候存在非常大的误差。
本质上是Euler 法只计算了每一步差值中的一阶部分，由Taylor级数可知：当公式中的二阶导数较大时就会产生明显的偏差，同时迭代过程中由于使用到上一部的结果，误差会在迭代中传播，因此这种Euler法在实际中是无法使用的，但是却给求解微分方程数值解提供了好的开始。
另外一种Euler法是隐式Euler法，其迭代公式是，它并没有解决上面所说的问题，同时它的计算更加繁琐，对于无法化简成显示迭代的公式时还需要用迭代法求解非线性方程。
为了解决上面的方法，就需要提高迭代公式中计算差值的阶数，下面介绍了梯形法和改进Euler法，它们都是二阶方法。
2.2.梯形公式和改进Euler法梯形公式以及改进Euler法都属于二阶方法，下面证明它是二阶方法，使用两次Taylor 公式，将和展开：将得到从上式可以看出，梯形法的局部截断误差的主要部分是，是关于步长的三次式，这说明了梯形法取到了差值中的二次项，因此梯形法是二阶方法。
从上面可以得到梯形的迭代公式：但是上式并不容易计算，因为上式中的为带求量，当无法化成显式形式时，需要对上式进行迭代求解。
因此梯形公式不易通过计算机编程求解，实际上改进的Euler法更容易求解。
改进Euler法迭代公式先通过显式Euler法求出一个估计值，通过这个估计值来计算，然后方程就变成了显式方程，从而可以得到修正值，改进Euler法更适合计算机编写程序，同样解决初值问题，详细MATLAB程序见附录2，得到的对比图像如图2.2所示。
图2.2 改进Euler法精确解与数值解对比由于改进Euler法用来求解的并不是精确解，所以得到的导数会有一定误差，因此改进Euler法的实际局部截断误差不仅仅是。
2.3.Euler法实用性从图2.1可以看出来一阶方法精确度非常差，基本上是无法用到实际工程中的，因此显式和隐式Euler法只是提供一种对微分方程求解的思想。
从图2.2中得到的数值解相对图2.1已经有了明显的改善，但是对于精度要求较高的工程问题，梯形法和改进Euler法这样的二阶方法同样是不满足要求的。
但通过这两个方法，了解到提高方法取差值中的更高阶数项能够达到提高精度的目的，后面内容中的Runge-Kutta法就是由此思想而来。
将细节部分放大更能够看出两种方法的精度，如图2.3所示(左图是Euler法，右图是改进Euler法)。
图2.3 两种方法细节部分3.Runge-Kutta Method(龙格库塔法)求解这一节的目的是弄清楚Runge-Kutta法(后面简称R-K法)的基本原理，并弄清楚MATLAB 中ode系列适合解决哪一种工程问题，希望达到的目的对于任何一个工程上提出的常微分方程问题都能够利用MATLAB求得误差在能够接受范围内的数值解。
3.1.Runge-Kutta基本原理MATLAB中使用R-K法的函数有ode23和ode45，其中ode是ordinary differential equation 的缩写，其中ode23表示使用的是2阶3级R-K法，ode45使用的是4阶5级R-K法。
将局部截断误差表示如下式：将上式中的方法称作阶级Runge-Kutta法。
其中：下面将构造出四阶四级经典R-K法，并通过MATLAB编程实现对初值问题求解。
四级四阶R-K法的本质是用四项级数来取得Taylor展开的前面四阶项，局部截断误差首项是步长的五次方，即：从上述表达式中可以看出，实际上R-K法是希望通过线性组合来近似Taylor公式中的阶导数组合，这样能够避免求解函数的阶导数，从而大量减小了计算量。
实际计算中将用来表示，则利用二元Taylor公式展开，舍去，并将同次数项系数相等，得到经典四阶四级R-K法的公式：上面从一阶到二阶的数值求法让误差有很大的改观，利用上述公式写出四阶四级的MATLABA程序(见附录3)，得到的图像如图3.1所示，图中解析解和数值解基本上重叠了，并且在处的误差已经小于了，而显式Euler法的误差为数量级，二阶改进Euler 法误差在数量级。
图3.1 四阶四级R-K法数值解和精确解对比3.2.MATLAB中使用Runge-Kutta法的函数MATLAB中利用R-K法求解常微分方程的函数主要是ode45，它用来求解非刚性常微分方程，是工程计算中首选函数。
另外的ode23是2阶3级R-K法函数，它可以用来求解轻微刚性方程，其精度相对较低，但是效率较高。
它们都是显式的R-K法，另外一个函数式ode23tb，它是隐式的R-K法可以用来求解非刚性的常微分方程。
4.使用MATLAB求解常微分方程MATLAB中求解常微分方程共有7个函数，表4-1列举出它们的使用范围以及简介。
表4-1 MATLAB中求解常微分方程函数通过表4-1了解到，当遇到一个常微分方程问题，首先使用ode45对其进行求解，如果求解非常慢或者无法求解，那么问题很可能是刚性问题，需要使用ode15s进行求解。
对于精度要求非常高的问题可使用ode113尝试进行求解，对于其他的刚性问题，需要具体研究使用哪一个函数。
4.1.使用ode45函数求解非刚性常微分方程ode45函数实现了变步长四阶五级Runge-Kutta-Felhberg算法(简称RKF法)。
调用格式如下所示。
公式中是函数的句柄，是求解区间的起点和终点，是初始值，后面的都是各种控制求解细节的选项。
其中的可以是匿名函数、子函数、嵌套函数、单独M 文件等形式。
实际上ode45能够求解常微分方程组，其中的使用状态变量表示方法，左侧是各个状态变量的导数，右侧是方程。
ode45可以解决高阶常微分方程的数值求解问题，基本思想是将高阶通过变量替换变成低阶，然后列成下面所示的微分方程组，然后进行求解。
对于处置初值问题，使用MATLAB求解(代码见附录4)如图3.2所示，其精度与上面自己编写的程序差不多，在处的误差小于。
图3.2 使用ode45求解细节4.2.刚性常微分方程如上所述，ode45函数能够求解的是非刚性问题，而刚性(stiff)问题使用ode45求解会产生求解缓慢，或者完全无法求解的情况，这种情况可以判断这个方程是刚性的。
可以考虑使用ode23或者ode23s函数进行求解。
刚性常微分方程指的是其特征根相差很大，但似乎没有非常有效判断一个方程是刚性还是非刚性的办法，需要实际求解的时候通过不同的函数进行求解比较。
5.总结文本主要研究了如何求常微分方程的数值解，从开始的一阶显式Euler法，到四阶R-K 法，精度提高了很多。
同时通过编写MATLAB代码来对具体初值问题进行求解，对课本上将的机制有了更加深入的了解，对具体算法过程有了深入的认识，从中学到了很多。
另外文中总结了MATLAB中求解常微分方程数值解函数的用法，在面对不同的工程问题的时候，怎样去挑选正确的函数来求解具体问题，更多的是学会了如何去解决实际问题。
实际问题中更多是由偏微分方程进行描述的，由于时间和个人数学水平原因，文中没有对偏微分方程数值求解进行研究，这是文章的不足，有待在后续的学习当中对其进行解决。
参考文献[1]于寅.高等工程数学[M].武汉:华中科技大学出版社.2012[2]李红.数值分析(第2版)[M].武汉:华中科技大学出版社.2010[3]胡建伟,汤怀民.微分方程数值方法(第二版)[M].北京:科学出版社.2007[4]薛定宇,陈阳泉.高等应用数学问题的MATLAB求解(第二版)[M].北京:清华大学出版社.2008附录1.显式Euler法数值求解clear,clc %y1中存储0~5中的数值解y2中存储-5~0中的数值解dyx=@(x,y) (-2*y-4*x); %方程右边函数fx=@(x) (exp(-2*x)-2*x+1); %解析解y1(1)=2;y2(1)=2;h = 0.1;i=2;for x=0:h:5-hy1(i)=y1(i-1)+h*dyx(x,y1(i-1));i =i +1;endi=2;for x=0:-h:-5+hy2(i) = y2(i-1)-h*dyx(x,y2(i-1));i=i+1;endx=-5:h:5;y=fx(x);ye=[fliplr(y2(2:end)) 2 y1(2:end)];plot(x,y,x,ye,'--');legend('\fontsize{12}\fontname{楷体}精确解',...'\fontsize{12}\fontname{楷体}数值解','interpreter','latex');text(x(8)+0.1,y(8),'$\leftarrow e^{-2x}-2x+1$','interpreter','latex','fontsize',15);2.改进Euler法数值求解clear,clc %y1中存储0~5中的数值解y2中存储-5~0中的数值解dyx=@(x,y) (-2*y-4*x); %方程右边函数fx=@(x) (exp(-2*x)-2*x+1); %解析解y1(1)=2;y2(1)=2;yy1(1)=2;yy2(1)=2; %yy1和yy2中存储的为估计值h = 0.1;i=2;for x=0:h:5-hyy1(i)=y1(i-1)+h*dyx(x,y1(i-1));yy1(i)=dyx(x+h,yy1(i));y1(i)=y1(i-1)+h/2*(dyx(x,y1(i-1))+yy1(i));i =i +1;endi=2;for x=0:-h:-5+hyy2(i) = y2(i-1)-h*dyx(x,y2(i-1));yy2(i) = dyx(x-h,yy2(i));y2(i) =y2(i-1)-h/2*(dyx(x,y2(i-1))+yy2(i));i=i+1;endx=-5:h:5;y=fx(x);ye=[fliplr(y2(2:end)) 2 y1(2:end)];plot(x,y,x,ye,'--');set(gcf,'color','white');legend('\fontsize{12}\fontname{楷体}精确解',...'\fontsize{12}\fontname{楷体}数值解','interpreter','latex');text(x(8)+0.1,y(8),'$\leftarrow e^{-2x}-2x+1$','interpreter','latex','fontsize',15);3.四阶四级Runge-Kutta法数值求解clear,clc %y1中存储0~5中的数值解y2中存储-5~0中的数值解dyx=@(x,y) (-2*y-4*x); %方程右边函数,即f(x,y)fx=@(x) (exp(-2*x)-2*x+1); %解析解y1(1)=2;y2(1)=2;h = 0.1;i=2;for x=0:h:5-hk1=dyx(x,y1(i-1));k2=dyx(x+h/2,y1(i-1)+h/2*k1);k3=dyx(x+h/2,y1(i-1)+h/2*k2);k4=dyx(x+h,y1(i-1)+h*k3);y1(i)=y1(i-1)+h/6*(k1+2*k2+2*k3+k4);i =i +1;endi=2;for x=0:-h:-5+hk1=dyx(x,y2(i-1));k2=dyx(x-h/2,y2(i-1)-h/2*k1);k3=dyx(x-h/2,y2(i-1)-h/2*k2);k4=dyx(x-h,y2(i-1)-h*k3);y2(i)=y2(i-1)-h/6*(k1+2*k2+2*k3+k4);i=i+1;endx=-5:h:5;y=fx(x);ye=[fliplr(y2(2:end)) 2 y1(2:end)];subplot(1,2,1);plot(x,y,x,ye,'--');set(gcf,'color','white');legend('\fontsize{12}\fontname{楷体}精确解',...'\fontsize{12}\fontname{楷体}数值解','interpreter','latex'); hold on;subplot(1,2,2);plot(x,y,x,ye,'--');legend('\fontsize{12}\fontname{楷体}精确解',...'\fontsize{12}\fontname{楷体}数值解','interpreter','latex'); axis([-5 -4.9995 2.2030e4 2.204e4]);4.使用ode45求解dyx=@(x,y) (-2*y-4*x); %方程右边函数fx=@(x) (exp(-2*x)-2*x+1); %解析解[xx,yy]=ode45(dyx,[0 5],2);[xxx,yyy]=ode45(dyx,[0 -5],2),x=-5:0.1:5;y=fx(x);sx=[flipud(xxx) xx];sy=[flipud(yyy) yy];subplot(1,2,1);plot(x,y,sx,sy,'--');subplot(1,2,2);plot(x,y,sx,sy,'--');set(gcf,'color','white');axis([-5 -4.9997 2.203e4 2.205e4]);。