下载文档原格式
微分方程的数值解法微分方程是自然科学和现代技术领域中一种最基本的数学描述工具，它可以描述物理世界中的各种现象。
微分方程的解析解往往很难求出，因此数值解法成为解决微分方程问题的主要手段之一。
本文将介绍几种常见的微分方程的数值解法。
一、欧拉法欧拉法是微分方程初值问题的最简单的数值方法之一，它是由欧拉提出的。
考虑一阶常微分方程：$y'=f(t,y),y(t_0)=y_0$其中，$f(t,y)$表示$y$对$t$的导数，则$y(t_{i+1})=y(t_i)+hf(t_i,y_i)$其中，$h$为步长，$t_i=t_0+ih$，$y_i$是$y(t_i)$的近似值。
欧拉法的精度较低，误差随着步长的增加而增大，因此不适用于求解精度要求较高的问题。
二、改进欧拉法改进欧拉法又称为Heun方法，它是由Heun提出的。
改进欧拉法是在欧拉法的基础上进行的改进，它在每个步长内提高求解精度。
改进欧拉法的步骤如下：1. 根据当前$t_i$和$y_i$估算$y_{i+1}$:$y^*=y_i+hf(t_i,y_i),t^*=t_i+h$2. 利用$y^*$和$t^*$估算$f(t^*,y^*)$:$f^*=f(t^*,y^*)$3. 利用$y_i$、$f(t_i,y_i)$和$f^*$估算$y_{i+1}$:$y_{i+1}=y_i+\frac{h}{2}(f(t_i,y_i)+f^*)$改进欧拉法具有比欧拉法更高的精度，但是相较于其他更高精度的数值方法，它的精度仍然较低。
三、龙格-库塔法龙格-库塔法是一种广泛使用的高精度数值方法，它不仅能够求解一阶和二阶常微分方程，还能够求解高阶常微分方程和偏微分方程。
其中，经典的四阶龙格-库塔法是最常用的数值方法之一。
四阶龙格-库塔法的步骤如下：1. 根据当前$t_i$和$y_i$估算$k_1$:$k_1=f(t_i,y_i)$2. 根据$k_1$和$y_i$估算$k_2$:$k_2=f(t_i+\frac{h}{2},y_i+\frac{h}{2}k_1)$3. 根据$k_2$和$y_i$估算$k_3$:$k_3=f(t_i+\frac{h}{2},y_i+\frac{h}{2}k_2)$4. 根据$k_3$和$y_i$估算$k_4$:$k_4=f(t_i+h,y_i+hk_3)$5. 根据$k_1$、$k_2$、$k_3$和$k_4$计算$y_{i+1}$:$y_{i+1}=y_i+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$龙格-库塔法的精度较高，在求解一些对精度要求较高的问题时，龙格-库塔法是一个比较好的选择。
微分方程的数值解法微分方程是描述自然界中众多现象和规律的重要数学工具。
然而，许多微分方程是很难或者无法直接求解的，因此需要使用数值解法来近似求解。
本文将介绍几种常见的微分方程数值解法。
1. 欧拉方法欧拉方法是最简单的数值解法之一。
它将微分方程转化为差分方程，通过计算离散点上的导数来逼近原方程的解。
欧拉方法的基本思想是利用当前点的导数值来估计下一个点的函数值。
具体步骤如下：首先，将自变量区间等分为一系列的小区间。
然后，根据微分方程的初始条件，在起始点确定初始函数值。
接下来，根据导数的定义，计算每个小区间上函数值的斜率。
最后，根据初始函数值和斜率，递推计算得到每个小区间上的函数值。
2. 龙格-库塔方法龙格-库塔方法是一种常用的高阶精度数值解法。
它通过进行多次逼近和修正来提高近似解的准确性。
相比于欧拉方法，龙格-库塔方法在同样的步长下可以获得更精确的解。
具体步骤如下：首先，确定在每个小区间上的步长。
然后，根据微分方程的初始条件，在起始点确定初始函数值。
接下来，根据当前点的导数值，使用权重系数计算多个中间点的函数值。
最后，根据所有中间点的函数值，计算出当前点的函数值。
3. 改进欧拉方法（改进的欧拉-克罗默法）改进欧拉方法是一种中阶精度数值解法，介于欧拉方法和龙格-库塔方法之间。
它通过使用两公式递推来提高精度，并减少计算量。
改进欧拉方法相对于欧拉方法而言，增加了一个估计项，从而减小了局部截断误差。
具体步骤如下：首先，确定在每个小区间上的步长。
然后，根据微分方程的初始条件，在起始点确定初始函数值。
接下来，利用欧拉方法计算出中间点的函数值。
最后，利用中间点的函数值和斜率，计算出当前点的函数值。
总结：微分方程的数值解法为我们研究和解决实际问题提供了有力的工具。
本文介绍了欧拉方法、龙格-库塔方法和改进欧拉方法这几种常见的数值解法。
选择合适的数值解法取决于微分方程的性质以及对解的精确性要求。
在实际应用中，我们应该根据具体情况选择最合适的数值解法，并注意控制步长以尽可能减小误差。
微分方程数值解使用数值方法求解微分方程微分方程是描述自然现象中变化的数学模型，它是数学和科学研究中的重要工具。
然而，许多微分方程并没有精确的解析解，因此需要使用数值方法来近似求解。
本文将介绍一些常用的数值方法来求解微分方程，包括欧拉方法、改进的欧拉方法和龙格-库塔方法。
一、欧拉方法欧拉方法是最简单、最基础的数值方法之一。
它基于微分方程解的定义，通过离散化自变量和因变量来逼近解析解。
假设我们要求解的微分方程为dy/dx = f(x, y)，初始条件为y(x0) = y0。
将自变量x分割成若干个小区间，步长为h，得到x0, x1, x2, ..., xn。
根据微分方程的定义，我们可以得到递推公式 yn+1 = yn + h*f(xn, yn)。
用代码表示即为：```def euler_method(f, x0, y0, h, n):x = [x0]y = [y0]for i in range(n):xn = x[i]yn = y[i]fn = f(xn, yn)xn1 = xn + hyn1 = yn + h*fnx.append(xn1)y.append(yn1)return x, y```二、改进的欧拉方法欧拉方法存在着局部截断误差，即在每个小区间上的误差。
改进的欧拉方法是对欧拉方法的改进，可以减小截断误差。
它的递推公式为yn+1 = yn + h*(f(xn, yn) + f(xn+1, yn+1))/2。
用代码表示即为：```def improved_euler_method(f, x0, y0, h, n):x = [x0]y = [y0]for i in range(n):xn = x[i]yn = y[i]fn = f(xn, yn)xn1 = xn + hyn1 = yn + h*(fn + f(xn1, yn + h*fn))/2x.append(xn1)y.append(yn1)return x, y```三、龙格-库塔方法龙格-库塔方法是一种更加精确的数值方法，它通过计算多个递推式的加权平均值来逼近解析解。
微分方程的数值解法微分方程（Differential Equation）是描述自然界中变化的现象的重要工具，具有广泛的应用范围。
对于一般的微分方程，往往很难找到解析解，这时候就需要使用数值解法来近似求解微分方程。
本文将介绍几种常见的微分方程数值解法及其原理。
一、欧拉方法（Euler's Method）欧拉方法是最基本也是最容易理解的数值解法之一。
它的基本思想是将微分方程转化为差分方程，通过给定的初始条件，在离散的点上逐步计算出函数的近似值。
对于一阶常微分方程dy/dx = f(x, y)，利用欧拉方法可以得到近似解：y_n+1 = y_n + h * f(x_n, y_n)其中，h是步长，x_n和y_n是已知点的坐标。
欧拉方法的优点在于简单易懂，但是由于是一阶方法，误差较大，对于复杂的微分方程可能不够准确。
二、改进的欧拉方法（Improved Euler's Method）改进的欧拉方法又称为改进的欧拉-柯西方法，是对欧拉方法的一种改进。
它通过在每一步计算中利用两个不同点的斜率来更准确地逼近函数的值。
对于一阶常微分方程dy/dx = f(x, y)，改进的欧拉方法的迭代公式为：y_n+1 = y_n + (h/2) * [f(x_n, y_n) + f(x_n+1, y_n + h * f(x_n, y_n))]相较于欧拉方法，改进的欧拉方法具有更高的精度，在同样的步长下得到的结果更接近真实解。
三、四阶龙格-库塔方法（Fourth-Order Runge-Kutta Method）四阶龙格-库塔方法是一种更高阶的数值解法，通过计算多个点的斜率进行加权平均，得到更为准确的解。
对于一阶常微分方程dy/dx = f(x, y)，四阶龙格-库塔方法的迭代公式为：k1 = h * f(x_n, y_n)k2 = h * f(x_n + h/2, y_n + k1/2)k3 = h * f(x_n + h/2, y_n + k2/2)k4 = h * f(x_n + h, y_n + k3)y_n+1 = y_n + (k1 + 2k2 + 2k3 + k4)/6四阶龙格-库塔方法是数值解法中精度最高的方法之一，它的计算复杂度较高，但是能够提供更为准确的结果。
微分方程是数学中的一种重要的方程类型，广泛应用于物理、工程、经济等领域。
解微分方程有各种方法，其中数值解法是一种重要而实用的方法。
微分方程的数值解法是通过数值计算来求解微分方程的近似解。
它的基本思想是将微分方程转化为差分方程，并用计算机进行迭代计算，从而求得微分方程的数值解。
数值解法的关键在于如何将微分方程转化为差分方程。
常见的方法有欧拉方法、改进欧拉方法、龙格－库塔方法等。
这些方法都是基于泰勒级数展开的原理进行推导的。
以欧拉方法为例，其基本思路是将微分方程中的导数用差商的方式近似表示，然后通过迭代计算，逐步逼近微分方程的解。
欧拉方法的具体步骤如下：首先确定微分方程的初始条件，即给定t0时刻的函数值y0，然后选取一定的步长ℎ，利用微分方程的导数计算差商y′=dy，进而根据差商dt得到下一个时刻的函数值y n+1=y n+ℎy′。
通过不断迭代计算，即可得到微分方程在一定时间区间内的数值解。
数值解法的另一个重要问题是误差控制。
由于数值计算本身的误差以及近似方法的误差，数值解法所得到的结果通常与真实解存在误差。
为了控制误差，常用的方法有缩小步长ℎ、提高近似方法的阶数等。
此外，还可以通过与解析解进行比较，评估数值解的准确性。
微分方程的数值解法具有以下几点优势。
首先，微分方程的解析解通常较难求得，而数值解法可以给出一个近似解，提供了一种有效的解决方案。
其次，数值解法可以利用计算机的高速运算能力，进行大规模复杂微分方程的求解。
此外，数值解法还可以在实际问题中进行仿真和优化，即通过调整参数来求解微分方程，从而得到最优解。
尽管微分方程的数值解法具有广泛的应用前景，但也存在一些问题和挑战。
首先，数值解法的稳定性和收敛性需要深入研究和分析。
其次，数值解法的计算量通常较大，对计算机运算能力和存储空间的要求较高。
此外，数值解法还需要对问题进行适当的离散化处理，从而可能引入一定的误差。
综上所述，“微分方程的数值解法”是一种重要而实用的方法，可以有效地求解微分方程的近似解。
微分方程数值解法微分方程是数学中的重要概念，它描述了物理系统中变量之间的关系。
解微分方程是许多科学领域中常见的问题，其中又可以分为解析解和数值解两种方法。
本文将重点介绍微分方程的数值解法，并详细讨论其中的常用方法和应用。
一、微分方程的数值解法概述微分方程的解析解往往较为复杂，难以直接求解。
在实际问题中，我们通常利用计算机进行数值计算，以获得方程的数值解。
数值解法的基本思想是将微分方程转化为一组离散的数值问题，通过逼近连续函数来获得数值解。
二、常见的数值解法1. 欧拉法欧拉法是最基础的数值解法之一，其核心思想是将微分方程转化为差分方程，通过逼近连续函数来获得数值解。
欧拉法的基本形式为：yn+1 = yn + h·f(xn, yn)其中，yn表示第n个时间步的数值解，h为时间步长，f为微分方程右端的函数。
欧拉法的精度较低，但计算简单，适用于初步估计或简单系统的求解。
2. 改进的欧拉法（Heun法）改进的欧拉法（Heun法）是对欧拉法的改进，其关键在于求解下一个时间步的近似值时，利用了两个斜率的平均值。
Heun法的基本形式为：yn+1 = yn + (h/2)·(k1 + k2)k1 = f(xn, yn)，k2 = f(xn+h, yn+h·k1)Heun法较欧拉法的精度更高，但计算量较大。
3. 龙格-库塔法（RK方法）龙格-库塔法是一类常用的数值解法，包含了多个不同阶数的方法。
其中，最常用的是经典四阶龙格-库塔法（RK4法），其基本形式为：k1 = f(xn, yn)k2 = f(xn + h/2, yn + (h/2)·k1)k3 = f(xn + h/2, yn + (h/2)·k2)k4 = f(xn + h, yn + h·k3)yn+1 = yn + (h/6)·(k1 + 2k2 + 2k3 + k4)RK4法实现较为复杂，但精度较高，适用于解决大多数常微分方程问题。
求微分方程数值解
微分方程数值解是一种数学方法，用于解决一些复杂的微分方程，特别是那些无法通过解析方法求解的微分方程。
通过数值解法，我们可以得到微分方程的近似解，并且可以在计算机上进行实现，以便更好地理解和分析问题。
我们需要将微分方程转化为差分方程，这样就可以利用数值方法进行求解。
差分方程是一种以离散形式表示微分方程的方法，通过近似替代微分表达式，将连续问题转化为离散问题，从而实现计算机求解。
常见的数值方法包括欧拉方法、龙格-库塔方法等，它们通过不断迭代求解差分方程，逼近微分方程的解。
在应用数值解法求解微分方程时，需要注意选择合适的步长和迭代次数，以确保数值解的准确性和稳定性。
步长过大会导致数值误差增大，步长过小则会增加计算量，影响计算效率。
因此，需要在准确性和效率之间寻找平衡点，选择合适的参数进行计算。
在使用数值解法时，还需要考虑边界条件和初值条件的设定。
这些条件对于微分方程的求解至关重要，不同的条件设定可能会导致不同的数值解，甚至无法得到有效的解。
因此，在进行数值计算之前，需要对问题进行充分的分析和理解，确定合适的条件，以确保数值解的准确性和可靠性。
总的来说，微分方程数值解是一种强大的工具，可以帮助我们解决
复杂的微分方程，探索未知的领域。
通过合理的数值方法和参数选择，我们可以得到准确的数值解，从而更好地理解和应用微分方程的理论。
希望通过不断的探索和实践，我们可以更深入地理解微分方程数值解的原理和方法，为科学研究和工程实践提供更多有益的帮助。
微分方程数值解差分法微分方程是自然科学和工程技术中广泛使用的工具，它们描述了许多物理过程的动力学行为。
对于复杂的微分方程，解析解往往很难或者不可能得到。
此时我们需要数值解差分法来解决问题。
一、微分方程数值解的方法1.分裂法分裂法是将一个复杂的微分方程分解为多个简单的方程。
例如，将一个偏微分方程分解成几个常微分方程，从而可以方便地使用数值方法计算解。
2.有限差分法有限差分法是一种常见的微分方程数值计算方法。
它将一维或多维的连续函数离散为一系列离散点，然后使用差分方程近似微分方程，最后用迭代法计算数值解。
3.有限元法有限元法是一种广泛使用的数值计算方法，它可以用于求解各种类型的微分方程。
该方法将求解区域分割成多个小区域，然后对每个小区域进行离散化和近似处理。
二、数值解差分法数值解差分法是微分方程数值解的基本方法之一。
它是一种基于差分方程的离散化方法，可以对微分方程进行近似，并将微分方程转化为一个差分方程。
数值解的差分法可以分为前向差分、后向差分和中心差分三种方法。
1.前向差分法前向差分法使用前一时间步的值，计算当前时间步的值。
它的近似误差随着时间步长的增大而增大。
前向差分的公式如下：y_i+1 = y_i + hf_i(x_i,y_i)其中，h是时间步长，f_i是微分方程的左侧。
2.后向差分法后向差分法使用后一时间步的值，计算当前时间步的值。
它的近似误差随着时间步长的增大而减小。
后向差分的公式如下：y_i+1=y_i + hf_i(x_i+1,y_i+1)3.中心差分法中心差分法使用前一时间步和后一时间步的值，计算当前时间步的值。
它的近似误差随着时间步长的增大而增大。
中心差分的公式如下：y_i+1=y_i + 1/2hf_i(x_i,y_i) + 1/2hf_i(x_i+1,y_i+1)三、差分法的优缺点差分法作为微分方程数值解的一种基本方法，具有以下优缺点：1.优点（1）简单易实现：差分法的实现很简单，只需要计算微分方程的离散值和靠近值即可。
微分方程的数值解法比较
微分方程是自然界中描述变化的数学工具之一。
在科学与工程领域, 微分方程的解对于预测系统行为、找到最优控制和优化过程等方面具有重要意义。
然而，通常情况下，微分方程很难通过解析的方式得到解析解。
因此，数值解法成为了解决微分方程的重要手段之一。
本文将介绍并比较几种常见的微分方程数值解法：
欧拉方法
欧拉方法是一种简单直观的数值解法，通过将一个微分方程转化为差分方程来近似求解。
其基本思想是通过微分方程在某一点的导数值来估计在该点邻近的函数值，然后逐步迭代以获得整体的近似解。
改进欧拉方法
改进欧拉方法在欧拉方法的基础上做出了改进，通过取相邻两点的斜率平均值来提高数值解的精度和稳定性。
龙格-库塔方法
龙格-库塔方法是一类广泛应用的数值解法，其基本思想是通过组合多个简单欧拉法来获得更高阶的精度，是当前应用最广泛的微分方程数值解法之一。
有限元方法
有限元方法是另一种常见的微分方程数值解法，主要用于处理偏微分方程。
通过将求解区域离散化为有限个单元，再进行有限维空间上的逼近来求解微分方程，其精度和收敛速度较高。
总结与比较
在实际应用中，不同的微分方程数值解法各有优缺点。
欧拉方法和改进欧拉方法简单易懂，但精度较低；龙格-库塔方法精度高但计算量大；有限元方法适用于处理复杂的偏微分方程，但较为复杂。
因此，在选择数值解法时需要根据具体情况综合考虑精度、计算量等因素。
以上就是一些常见微分方程数值解法的介绍，希望能对读者有所帮助。
微分方程的常用数值解法摘要：微分方程是数学中的一种重要的方程类型，它能描述自然现象和工程问题中的许多变化规律。
但是大多数微分方程解法是无法用解析的方式求解的，因此需要借助数值解法来近似求解。
本文将介绍微分方程的常用数值解法。
关键词：欧拉方法；龙格-库塔方法；微分方程；常用数值解法一、微分方程数值解方法微分方程数值解法是数学中的重要部分。
欧拉方法、龙格-库塔方法和二阶龙格-库塔方法是常用的微分方程数值解法，下面就分别介绍这三种方法。
（一）欧拉方法欧拉方法是解初值问题的一种简单方法，它是欧拉用的第一种数值方法，也叫向前欧拉法。
欧拉方法是利用微分方程的定义式y’=f(x, y)，将它带入微分方程初值问题y（x_0）=y_0中，以y_0为初始解，在每一步上通过沿着切线的方法进行估计并推进新的解y_{i+1}：y_i+1=y_i+hf(x_i,y_i)其中，x_i和y_i是我们知道的初始条件，h是求解过程中的步长，f是微分方程右端项。
它是一种时间迭代的算法，易于实现，但存在着精度不高的缺点。
（二）龙格-库塔方法龙格-库塔方法是一种经典迭代方法，也是近代微分方程数值解法发展的里程碑之一。
龙格-库塔方法的主要思想是利用规定的阶码及阶向量，通过递推求解微分方程数值解的近似值。
龙格-库塔方法的方式不同，其步骤如下：第一步：根据微分方程，计算出在x_i和y_i的值。
第二步：在x_i处对斜率进行估计，并利用这个斜率来求解下一步所需的y_i+1值。
第三步：使用x_i和y_i+1的值来重新估计斜率。
第四步：使用这个新的斜率来更新y_i+1的值。
（三）二阶龙格-库塔方法二阶龙格-库塔方法是龙格-库塔方法的一种变体，它根据龙格-库塔方法的思想，使用更好的步长来提高数值解的精度。
二阶龙格-库塔方法的基本思路是，在第一次迭代时使用一个阶段小一半的y_i+1，然后使用这个估算值来计算接下来的斜率。
通过这种方法，可以提高解的精度。
二阶龙格-库塔方法的步骤如下：第一步：计算出初始阶段的y_i+1值。
微分方程的数值解与数值方法微分方程是数学中的重要内容，它描述了许多自然现象和物理问题中的变化规律。
解微分方程是求解已知条件下未知函数的问题，是数学建模和科学研究中的核心内容之一。
传统的解微分方程的方法有解析解和数值解两种，解析解是通过推导和运算得到的精确解，而数值解是通过近似计算获得的近似解。
本文将介绍微分方程的数值解方法和数值解的优缺点。
微分方程的数值解方法主要有两种：欧拉方法和改进的欧拉方法。
欧拉方法是一种基本的数值解方法，它根据微分方程在某一点的斜率来近似计算下一个点的函数值。
具体来说，欧拉方法将微分方程中的导数用差商表示，然后根据差商计算下一个点的函数值。
欧拉方法的优点是简单易懂，容易实现。
缺点是精度较低，容易产生误差。
改进的欧拉方法是对欧拉方法的改进，它通过考虑两个相邻点的斜率的平均值来计算下一个点的函数值。
改进的欧拉方法相对于欧拉方法来说，精度更高，误差更小。
数值解的优点是能够得到近似解，可以在一定程度上对实际问题进行模拟和仿真。
数值解方法对于复杂的微分方程或者无法求得解析解的微分方程非常有用。
数值解还可以帮助研究者验证解析解的正确性，并且可以用于求解一些实际问题，如物理问题和工程问题。
数值解的缺点是精度不如解析解高，容易产生误差，并且对初始条件和步长敏感。
此外，数值解的计算量较大，需要使用计算机来实现，而解析解则可以通过手工计算得到。
数值解方法在实际应用中有广泛的应用。
例如，微分方程在物理学中的应用非常广泛，如运动学和力学中的运动方程、电磁学中的麦克斯韦方程、量子力学中的薛定谔方程等。
这些方程往往是复杂的，无法通过解析方法求得精确解，只能通过数值解方法进行求解。
另外，数值解方法也在生物学、经济学、地理学等领域有重要的应用。
生物学中的生物动力学方程、经济学中的经济增长方程、地理学中的模拟气候变化等问题都需要通过数值解方法求解。
总结起来，微分方程的数值解方法是一种求解微分方程的有效工具。
第四章 微分方程数值解4.1 算 法当常微分方程能解析求解时，可利用Matlab 符号工具箱中的功能找到精确解. 见下例求解方程20y y y '''++=. 键入：syms x y %定义符号变量 diff_equ= ‘D2y+2*Dy -y=0’; %D2y 表示,y ''Dy=y 'y=dsolve (diff_equ, ‘x’)%定义x 为自变量y=cl*exp ((2^(1/2)-1)*x+c2*exp (-(2^(1/2)+1)*x) %表达式中含c1与c2，表示通解.%初始条件为y (0)=0,y '(0)=1时，按如下方式调用y=dsolve (diff_equ, ‘y (0)=0’, ‘Dy (0)=1’, ‘x’)y=1/4*2^(1/2)*exp ((2^(1/2)-1)*x)— 1/4*2^(1/2)*exp (-(2^(1/2)+1)*x) %画出函数y=y (x)的图形ezplot (y,[-2,2])图形具体形式请上机试之.在方程无法获得解析解的情况下，可方便地获得数值解. 下面的例子说明用Matlab 求数值解的方法及应注意的问题.[例1] 求解范德堡（vander pol ）方程222(1)0d x dx x x dt dtμ--+=求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，实现这一变换1,2/y x y dx dt ==则令 1/2dy dt y =22/(11)*21dy dt y y y μ=--编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为<待求解方程的函数名.m >，第二部分为求解主程序，本例中取名为main1.m.首先编写待求解方程的M -文件. 文件存盘名为“vdpol.m ”. function yprime=vdpol (,)t y yprime (1)=y (2);mu=2yprime (2)=mu*(1(1)^2)*(2)(1)y y y --;yprime=[yprime (1);yprime (2)];说明 函数yprime=vdpol (,)t y 中. 定义t 为自变量，y 的形式取决于求解方程的阶数，本例中，[(1),(2)],(1)y y y y =为解向量，(2)y 为导数向量. yprime (1)(1)y '=， yprime (2)(1)y ''=，函数返回vander pol 方程的导数列向量. 因为所求结果为方程数值解，所以各向量维数只有在主程序求解时定下精度后才能确定.主程序定名为main1.m ，你可用你所喜欢的其它名子，但vdpol.m 除外. clear functions%调试程序时，放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序 [,t y ]=ode23 (‘vdpol’,[0,30],[1,0]); y1=y (:,1); %解曲线. y2=y (:,2); %解曲线的导数. polt (,1,,2,t y t y ‘_ _’)说明 龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为： [,t x ]=ode23 (‘f ’,,0,ts x options) [,t x ]=ode45 (‘f ’,,0,ts x options)其中ts 可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. （1）若令[0,1,,]ts t t tf =，则输出在指定时刻0,1,,t t tf 给出，当0::ts t k tf =时，输出在区间[0,]t tf 的等分点上给出，k 为步长.（2）若[0,],0ts t tf t =为自变量初值，tf 为终值，此时，options 决定自变量t 的维数，t 中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长. 用于设定误差限的参数options 可缺省，此时系统设定相对误差为310-，绝对误差为610-，若自行设定误差限，可用如下语句： options=odeset (‘reltol’,rt , ‘abstol’,at ) 这里的rt 与at 分别为设定的相对与绝对误差.须注意的是无论用哪种方法确定t 的取值方式，0,t tf 必须由使用者确定且应与0x 相匹配. 0x 为初始条件，本例中0[1,0]y =，因为[0,30]ts =，这意味着解曲线(0)1y =，(0)0y '=，一般说，当解n 个未知函数的方程组时，0x 为n 维向量，共含有n 个初始条件.两个输出参数是列向量t 与矩阵x ，它们具有相同的行数，而矩阵x 的列数等于方程组的个数，本例中y 的列数为2，其中，(:,1)y 为自变量t 上各点函数值，(:,2)y 为t 上各点导数值.最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法.4.2 欧拉与龙格-库塔方法设有一阶方程与初始条件0(,)()y f x y y x y '=⎧⎨=⎩ （4.1） 其中(,)f x y 适当光滑，关于y 满足Lipschitz 条件，即存在L 使1212(,)(,)f x y f x y L y y -≤-则（4.1）式的解存在且惟一.关于()y y x =的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点12n x x x <<<<上寻求其数值近似解12,,,,n y y y .相邻两个节点间的间距1n n n h x x +=-称为步长，一般地取等步长h ，则0,1,2,n x x nh n =+=1、欧拉方法在区间1[,]n n x x +上用差商1()()n n y x y x h+-代替（4.1）式中y '，对(,)f x y 中x 在1[,]n n x x +上取值n x 还是1n x +，而形成向前欧拉公式与向后欧拉公式.（1）向前欧拉公式(,)f x y 取左端点n x ，得如下公式1()()(,())n n n n y x y x hf x y x +≈+ （4.2）从0x 点出发，由初值00()y x y =代入（4.2）求得1000(,)y y hf x y =+ （4.3）反复利用（4.2），有1(,) 0,1,2,n n n n y y hf x y n +=+= （4.4）其几何意义如图4.1所示.图中()y y x =为方程（4.1）的精确 解曲线，其上任意点(,)x y 处切线斜率为(,)f x y . 从初值000(,)P x y 点出发，用该点斜率00(,)f x y 作一直线段，在1x x = 处得到111(,)P x y ，1y 由（4.2）式确定， 再从111(,)P x y 出发，以11(,)f x y 为斜率 作直线段，在2x x =处得到222(,)P x y ，y0yO0x 误差4x1x 2x 3x x()y y x =3()y x32()y y x - 0P1P2P3P 4P3y这一过程继续下去，形成折线012P PP ，作为积分曲线()y y x =的近似，用1()n y x + 图4.1 表示在1n x +处的精确值，1n y +为解的近似值，不难得到23211()()()()2n n n h y x y y x O h O h ++''-=+=这一误差称为局部截断误差. 若一种算法局部截断误差为1()P O h +，则称该算法具有P 阶精度，所以向前欧拉公式具有1阶精度.（2）向后欧拉公式若(,)f x y 中x 取1[,]n n x x +中的1n x +，则有如下公式：111(,) 0,1,2,n n n n y y hf x y n +++=+= （4.5）称式（4.5）为向后欧拉公式，因为此式中1n y +未知，故称其为隐式公式，无法用其直接计算1n y +，一般用向前欧拉公式产生初值.(0)11(,) 0,1,2,n n n n y y hf x y n ++=+=再按下式迭代(1)()111(,),0,1,,0,1,k k n n n n y y hf x y k n ++++=+==其误差估计如下23211()()()()2n n n h y x y y x o h o h ++''-=+=精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2. 为讨论局部截断误差，在图4.2中设点 (,)n n n P x y 落在积分曲线()y y x =上，按式 （4.4）及式（4.5）分别得1n P +点为A 与B ， 且,A B 点一定在积分曲线()y y x =上相应 点θ的上、下两边,所以将式（4.4）与（4.5） 平均之，一定能得到更好的结果. （3）梯形公式 将向前与向后欧拉公式加以平均得到所 图4.2谓梯形公式111[(,)(,)] 0,1,2,2n n n n n n hy y f x y f x y n +++=++= （4.6）其局部截断误差为3()O h ，具有2阶精度.（4）改进的欧拉公式为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步1(,)n n n n y y hf x y +=+x 1n x +n x y()y y x =n PAθB111[(,)(,)], 0,1,2,2n n n n n n hy y f x y f x y n +++=++= （4.7）或写为1121211()2(,)(,)n n n n n n h y y k k k f x y k f x y hk ++⎧=++⎪⎪=⎨⎪=+⎪⎩0,1,2,n= （4.8）最后指出，上述欧拉方法可推广至微分方程组，如0000(,,)(,,)()()y f x y z z g x y z y x y z x z '=⎧⎪'=⎪⎨=⎪⎪=⎩ 向前欧拉公式为11(,,)(,,)n n n n n n n n n n y y hf x y z z z hg x y z ++=+⎧⎨=+⎩ 0,1,2,n=2、龙格_库塔方法 由微分中值定理1[()()]/(),01n n n y x y x h y x h θθ+'-=+<<又因为(,)y f x y '=，所以()(,())n n n y x h f x h y x h θθθ'+=++从而有1()()(),()n n n n y x y x hf x h y x h θθ+=+++ （4.9）令(,())n n k f x h y x h θθ=++，称其为区间1[,]n n x x +上的平均斜率，由（4.9）可知，给出一种平均斜率，可相应导出一种算法. 向前欧拉公式中(,)n n k f x y =，精度低. 改进欧拉公式中取111((,)(,))2n n n n k f x y f x y ++=+，精度提高，下面，我们在区间1[,]n n x x +内多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果. （1）2阶龙格_库塔公式11122121()(,)(,),0,1n n n n n n y y h k k k f x y k f x ah y hk a λλββ+=++⎧⎪=⎨⎪=++<<⎩ （4.10） 其中12211,,12a aβλλλ+===，由于4个未知数只有3个方程，所以解不惟一，若令121,12a λλβ+===，即得改进的欧拉公式，具有2阶精度.（3）4阶龙格_库塔公式 只给出精典格式中最常用的一种.112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y k k k k k f x y h h k f x y k h h k f x y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （4.11） 其计算精度为4阶4.3 模型与实验1、模型与问题 [例2] 单摆运动图4.3中一根长l 的细线，一端固定，另一端悬挂质量为m 的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度θ，然后使其自由运动，在不考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动. 0θ=为平衡位置，在小球摆动过程中，当与平衡位置夹 角为θ时，小球所受重力在其动运轨迹的分量为sin mg θ- （负号表示力的方向使θ减少），由牛顿第二定律可得微分方 程()sin ml t mg θθ''=- （4.12） 设小球初始偏离角度为0θ，且初速为0，式（4.12）的初 始条件为0(0),(0)0θθθ'== （4.13）当0θ不大时，sin θθ≈，式（4.12）化为线性常系数微分方程 图4.30g lθθ''+= （4.14）解得θlθmg0()t θθ= （4.15）简谐运动的周期为2T =. 现在的问题是：当0θ较大时，仍用θ近似sin θ，误差太大，式（4.12）又无解析解，试用数值方法在0030,10θθ==两种情况下求解，画出()t θ的图形，与近似解（4.15）比较，这里设25cm l =. [例3] 捕食与被捕食当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻t ，小鱼的数量为()x t ，鲨鱼的数量为()y t ，当甲独立生存时它的（相对）增长率与种群数量成正比，即有()()x t rx t '=，r 为增长率，而乙的存在使甲的增长率r 减少，设减少率与乙的数量成正比，而得微分方程()()(())x t x t r ay t rx axy '=-=- （4.16）比例系数a 反映捕食者掠取食饵的能力.乙离开甲无法生存，设乙独自存在时死亡率为d ，()()y t dy t '=-，甲为乙提供食物，使乙的死亡率d 降低，而促其数量增长，这一作用与甲的数量成正比，于是()y t 满足()(())y t y d bx t dy bxy '=-+=-+ （4.17）比例系数b 反映甲对乙的供养能力，设若甲，乙的初始数量分别为00(0);(0)x x y y == （4.18） 则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量()x t 、()y t 随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题：（1）设001,0.5,0.1,0.02,25,2r d a b x y ======，求方程（4.16），（4.17）在条件（4.18）下的数值解，画出(),()x t y t 的图形及相图(,)x y ，观察解(),()x t y t 的周期变化，近似确定解的周期和,x y 的最大、小值，近似计算,x y 在一个周期内的平均值. （2）从式（4.16）和（4.17）消去dt 得到()()dx x r ay dy y d bx -=-+ （4.19） 解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数. （3）将方程（4.17）改写为1()[]y x t d b y'=+ （4.20）在一个周期内积分，得到()x t 一周期内的平均值，类似可得()y t 一周期内的平均值，将近似计算的结果与理论值比较. 2、实验（1）方程（单摆问题）0sin (0),(0)0ml mg θθθθθ''=-⎧⎨'==⎩无解析解，为求其数值解，先将其化为等价的一阶方程组. 令12,x x θθ'==，方程化为1221102sin (0),(0)0x x g x x l x x θ'=⎧⎪⎪'=-⎨⎪==⎪⎩ 其中，09.8,25,g l θ==为100.1745≈（弧度）与300.5236≈（弧度）两种情况，具体编程如下：先以danbai.m 为文件名存放待解方程. 键入： function xdot =danbai (t,x) %x=[x (1),x (2)],g=9.8;1=25; %x (1)为解向量, x (2)是导数 xdot (1)=x (2); %xdot (1)=x '(1)=θ' xdot (2)=-g/1*sin (x (1)); %xdot (2)=θ''xdot=[xdot (1);xdot (2)]; %必须将导数向量以列向量形式给出.再以主程序(求数值解)调用待求方程，主程序用main2.m 为文件名存盘，其代码如下： clear functions[t,x]=ode23 (‘danbai’,[0,10],[0.1745,0])%只计算0100.1745θ==（弧度）的情形. %对近似解0()cos ,t wt w θθ==2T = w=sqrt (9.8/25); y=0.1745*cos (w*t);[t,x (:,1),y] %显示数据，无分号. hold on %欲在同一幅图中画近似解. plot (t,x (:,1), ‘b’) %画数值解，绿色 plot (t,y, ‘g*’) %用*号，红色画近似解. Hold off（2）食饵_捕食者 方程()x t rx axy '=- ()y t dy bxy '=-+可化为如下形式00x r ay x d bx y y ⎡⎤-⎡⎤⎡⎤⎢⎥=⎢⎥⎢⎥⎢⎥-+⎣⎦⎣⎦⎣⎦初始条件00(0),(0)x x y y ==表示为00(0)(0)x x y y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦以shier.m 存盘如下代码 function xdot=shier (t,x) r=1;d=0.5;a=0.1;b=0.02;xdot=diag ([r-a*x (2),-d+b*x (1)])*x;%x=(1)(2)x x x y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦,xdot=x y ⎡⎤⎢⎥⎢⎥⎣⎦以main3.m 存盘如下代码. clear functions ts=0:0.1:15; x0=[25,2];[t,x]=ode45 (‘shier’,ts,x0);[t,x] %显示数据t,x,y plot (t,x)gird %加网格线 gtext (‘x(t)’),gtext (‘y(t)’), %用点鼠标方式pause,figure (2) %将x 1(t),x 2(t)放至指定点 plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图xlabel (‘x’),ylabel (‘y’)可以猜测，1()x t 与2()x t 是周期函数，相图12(,)x x 是封闭曲线，从数值解可近似定出周期为110.7,x 的最大和最小值分别为99.3与22.0,x 的最大和最小值分别为28.4和2.0，为求1x 与2x 在一个周期的平均值，只需键入： y1=x (1:108,1); %x 1周期为10.7. x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积 x2p=trapz (y2)*0.1/10.7, %分值10711()1(1)2i i y i y i x =++∆∑ 可得1225,10x x == 对方程()()dx x r ay dy y d bx -=-+化为d bx r aydx dy x y-+-=积分得解()()d bx r ay x e y e c --= （4.21）即为原方程组的相轨线，其中c 由初始条件确定. 为说明上述相轨线是封闭的，令();()d bx r ay f x x e g y y e --==设其最大值分别为,m m f g ，若00,x y 满足00(),()m m f x f f y g ==则有00,d rx y b a==（令0,0f g ''==，可解出00,x y ），又当m m f g c ⋅≥时，相轨线（4.21）有意义. 当m m f g c =时，相轨线退化为一个点00(,)P x y ，对0m m c f g <<时，相轨线如图 4.4（c ），而图（a ），（b ）分别为()f x 与()g y 的图像，下面讨论如何由（a ），（b ）画（c ）.（a ） （b ） （c ）图4.4设(0)m m c pg p f =<<，若令0y y =，则有()f x p =，由（a ）知，12,x x ∃，使12()()f x f x p ==，且102x x x <<于是相轨线应过110(,)q x y 与220(,)q x y ，对12(,)x x x ∈，()f x p >，由()()()m m f x g x pg g y g =⇒<，令()g y q =，又由（b ）知，存在12,y y 使12()()g y g y q ==，且102y y y <<，于是相轨线又过31(,)q x y 与42(,)q x y 两点，所以对12,q q 间每个x ，轨线总要通过纵坐标为12102,()y y y y y <<的两点，于是相轨线是一条封闭曲线.相轨线封闭等价于(),()x t y t 是周期函数，记周期为T ，为求其在一个周期内的平均值(,)x y ，由1()()yx t d b y=+两边在一个周期内积分有((0)())y y T =：011ln ()ln (0)()T y T y dT d x x t dt T T b b b-⎡⎤==+=⎢⎥⎣⎦⎰ 同理()f xmfP1x 0x 2xxyy 1y ()()n m n mf x fg y g ==2yyqm g()g vx1y 0y 2y1xx 2x1q 3q2q00(,)P x yr y a=从而 00,x x y y ==即(),()x t y t 的平均值恰为相轨线中心点p 的坐标，提请读者注意的是：这里的0x 与0y 与初始条件中的00,x y 不是一件事.注意到,,,r d a b 在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者的能力成反比.3、练习内容（1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较.a ）22,(0)0y x y y '=-=或(0)1y =.b)222()0x y xy x n y '''++-=2()2,()22y y πππ'==-（Bessel 方程，这里令12n =，其精确解为sin y =. c)cos 0,(0)1,(0)0y y x y y '''+===. （2）倒圆锥形容器，上底面直径为1.2m ，容器的高亦为1.2m ，在锥尖的地方开有一直径为3cm 的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为h 时，水从小孔中流出的速度为v g =为重力加速度，若孔口收缩系数为0.6（即若一个面积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？（3）一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸上B 点，已知河水流速1v 与船在静水中的速度2v 之比为k .（a ）建立小船航线的方程，求其解析解.（b ）设12100m,1m/s,2m/s d v v ===，用数值解法求渡河所需时间，任意时刻小船的位置及航行曲线，作图并与解析解比较.（c ）若流速1v 为0,0.5,2 (m/s)，结果将如何？（4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律1212()(1),()(1)x y x t r x y t r y n n =-=- 其中，(),()x t y t 分别为t 时刻甲，乙两个种群的数量，12,r r 为其固有增长率，12,n n 为它们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长产生影响，将甲的方程修改为1112(1)x y x r x s n n =-- （4.22） 这里1s 的含意是：对于供养甲的资源而言，单位数量乙（相对2n ）的消耗率为单位数量甲（相对1n ）消耗的1s 倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为2212()(1)x y y t r y s n n =-- （4.23） 给定种群的初始值为 00(0),(0)x x y y == （4.24）及参数121212,,,,,r r s s n n 后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析解不存在，试用数值解法研究以下问题：（a ）设121212001,100,0.5,2,10r r n n s s x y ========，计算(),()x t y t ，画出它们的图形及相图(,)x y ，说明时间t 充分大以后(),()x t y t 的变化趋势（人们今天看到的已经是自然界长期演变的结果）.（b ）改变1200,,,r r x y ，但1s 与2s 不变（保持121,1s s <>），分析所得结果，若121.5(1),0.7(1)s s =>=<，再分析结果，由此你得到什么结论，请用各参数生态学上的含义作出解释.（c ）试验当120.8(1),0.7(1)s s =<=<时会有什么结果；当121.5(1), 1.7(1)s s =>=>时又会出现什么结果，能解释这些结果吗？。
微分方程的数值解法微分方程是数学中的一种重要的基础理论，广泛用于科学技术的研究中。
微分方程的解析解往往比较难求得，而数值解法则成为了解决微分方程的重要手段之一。
本文将阐述微分方程的数值解法，探讨一些经典的数值方法及其应用。
一、数值解法的基本思想微分方程的数值解法的基本思想是建立微分方程的差分方程，然后通过数值计算的方法求得差分方程的近似解，最终得到微分方程的数值解。
其中，差分方程是微分方程的离散化，将微分方程转化为差分方程的过程称为离散化或网格化。
离散化的目的是将连续问题转化为离散问题，使问题求解更为方便。
差分方程的计算通常需要将区间分成若干份，每一份都对应着一个节点，节点的个数与区间长度有关。
在每个节点处采集函数值，根据这些函数值计算出差分方程的值，再根据差分方程的迭代公式计算出每个节点的函数值。
因此差分方程的求解问题就转化成了求解节点函数值的问题。
二、欧拉法欧拉法是微分方程数值解法中最简单的一种方法，广泛应用于各种领域。
欧拉法的基本思想是运用泰勒公式，将函数在某一点展开成一次多项式，用两个相邻节点之间的差分来逼近导数的值，从而得到连续问题的离散解。
具体实现过程如下：1. 将微分方程的初始值问题区间[a,a]分成若干个小区间，每个小区间长度为a，共有a个节点，其中节点序列为a0,a1,a2,⋯,aa，节点之间的间隔为a。
2. 根据微分方程的迭代公式得到差分方程，即令aa+1=aa+aa(aa,aa)3. 按照差分方程的迭代公式，从初始值a0开始，逐一计算得到函数值，a1,a2,⋯,aa。
欧拉法的精度比较低，误差常常会较大，但是它运算速度快，实现简单，计算量小，因此在计算简单模型时常常使用。
三、龙格-库塔法龙格-库塔法是微分方程数值解法中精度最高的一种方法，具有比欧拉法更精确、更稳定的特点，广泛应用于各种实际问题中。
龙格-库塔法的主要思想是用多阶段逼近法估算每一步的函数值，从而提高时间的精度。
具体实现过程如下：1. 将微分方程的初始值问题区间[a,a]分成若干个小区间，每个小区间长度为a，共有a个节点，其中节点序列为a0,a1,a2,⋯,aa，节点之间的间隔为a。
微分方程数值解法概述微分方程是描述自然界和社会科学中许多现象的重要数学模型，它们在科学研究和工程技术中具有广泛的应用。
为了求解微分方程，人们开发了多种数值解法。
本文将对微分方程数值解法进行概述，介绍其中常用的几种方法。
一、欧拉方法欧拉方法是最简单的一种数值解法，它基于微分方程的定义。
欧拉方法将微分方程的解曲线离散化为一系列连接相邻点的线段，并通过计算斜率来近似曲线的切线。
具体步骤如下：1. 将解曲线上的点等距离地选取为x0, x1, x2, ..., xn。
2. 根据微分方程得出差分方程：y_(k+1) = y_k + h * f(x_k, y_k)，其中h为步长。
3. 通过迭代计算，得到近似解的数值解。
尽管欧拉方法简单直观，但由于是一阶方法，它的精度相对较低，容易出现截断误差。
二、改进的欧拉方法为了提高数值解的精度，人们改进了欧拉方法，例如改进的欧拉方法、改进的欧拉法和四阶改进的欧拉法等。
这些方法主要通过引入更高阶的项来减小截断误差，从而提高数值解的精度。
其中最常用的是四阶改进的欧拉法，也称为四阶龙格-库塔法（RK4）。
该方法具体步骤如下：1. 根据微分方程的定义，设置初始值y0。
2. 根据微分方程，计算中间点的斜率k1，k2，k3和k4。
3. 计算步长h * (k1+2k2+2k3+k4)/6，得到下一个节点的近似解。
4. 重复步骤2和步骤3，直到得到满足要求的数值解。
三、龙格-库塔方法（RK方法）龙格-库塔方法是一类经典的数值解法，常用于求解常微分方程。
与欧拉方法和改进的欧拉方法不同，龙格-库塔方法中的每个节点都有自己的权重。
最常用的是四阶龙格-库塔方法（RK4），其步骤与上述四阶改进的欧拉法类似。
四、有限差分法有限差分法是求解微分方程的一种常见数值方法。
该方法将微分方程中的导数用差商的形式进行近似，然后通过在离散的网格点上求解代数方程组来得到数值解。
有限差分法的核心思想是使用差商来逼近导数。
微分方程数值解法微分方程数值解法是一种将微分方程的解转化为数值计算的方法。
常用的微分方程数值解法包括欧拉法、隐式欧拉法、龙格-库塔法等。
1. 欧拉法：欧拉法是最简单的一种数值解法，它基于微分方程的定义，在给定的初始条件下，通过不断迭代计算微分方程在给定区间上的近似解。
欧拉法的迭代公式为：y_{n+1}=y_n+h\\cdot f(t_n,y_n)，其中y_n表示第n步的近似解，t_n表示第n步的时间，h表示步长，f(t_n,y_n)表示微分方程的右侧函数。
2. 隐式欧拉法：隐式欧拉法是欧拉法的改进，它在计算近似解时使用了未知公式的近似值，从而提高了精度。
隐式欧拉法的迭代公式为：y_{n+1}=y_n+h\\cdotf(t_{n+1},y_{n+1})，其中y_{n+1}表示第n+1步的近似解，t_{n+1}表示第n+1步的时间，h表示步长，f(t_{n+1},y_{n+1})表示微分方程的右侧函数。
3. 龙格-库塔法：龙格-库塔法是一种常用的高阶数值解法，它通过计算微分方程的斜率来提高精度。
最常见的是四阶龙格-库塔法，它的迭代公式为：y_{n+1}=y_n+\\frac{1}{6}(k_1+2k_2+2k_3+k_4)，其中k_1=h\\cdot f(t_n,y_n)，k_2=h\\cdotf(t_n+\\frac{h}{2},y_n+\\frac{1}{2}k_1)，k_3=h\\cdotf(t_n+\\frac{h}{2},y_n+\\frac{1}{2}k_2)，k_4=h\\cdotf(t_n+h,y_n+k_3)。
这些方法的选择取决于问题的性质和精度要求。
其中，欧拉法是最简单的方法，但精度较低，龙格-库塔法精度较高，但计算量较大。
在实际应用中需要根据问题的具体情况选择合适的数值解法。
微分方程数值解及其应用绪论 自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和科学技术中的实际问题的研究中, 常常需要求解微分方程．但往往只有少数较简单和典型的微分方程可求出其解析解，在大多数情况下,只能用近似法求解，数值解法是一类重要的近似方法．本文主要讨论一阶常微分方程的初值问题的数值解法，探讨这些算法在处理来自生活实际问题中的应用，并结合MATLAB 软件，动手编程予以解决． １ 微分方程的初值问题[1]1.1 预备知识在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题00(,)()dy f x y dx y x y ⎧=⎪⎨⎪=⎩ （1）这里(),f x y 是矩形区域R ：00,x x a y y b -≤-≤上的连续函数．对初值问题（1）需要考虑以下问题：方程是否一定有解呢?若有解，有多少个解呢?下面给出相关的概念与定理．定义1 Lipschitz 条件[1][2]：矩形区域R ：00,x x a y y b -≤-≤上的连续函数(),f x y 若满足：存在常数0L >，使得不等式()()1212,,f x y f x y L y y -≤-对所有()()12,,,x y x y R ∈都成立，则称(),f x y 在R 上关于y 满足Lipschitz 条件. 定理 1 解的存在唯一性定理[1][3]：设f 在区域()}{,,D x y a x b y R =≤≤∈上连续，关于y 满足Lipschitz 条件，则对任意的[]00,,∈∈x a b y R ，常微分方程初值问题（1）当[],x a b ∈时存在唯一的连续解()y x ．该定理保证若一个函数(),f x y 关于y 满足Lipschitz 条件，它所对应的微分方程的初值问题就有唯一解.在解的存在唯一性得到保证的前提下，自然要考虑方程的求解问题．求解微分方程虽然有多种解析方法，但根据工程和科学实践问题所得到的微分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，也往往因形式过于复杂或计算量太大而不实用，因此从实际问题中归结出来的微分方程主要依靠数值解法．定义 2 微分方程数值解：对初值问题（1）寻求数值解就是寻求解()y x 在一系列离散节点上的近似解0121,,,,,,n n y y y y y +L L ，相邻两个节点的间距1n n n h x x +=-称为步长.在一般情况下假定()0,1,i h h i ==L 为常数，这时节点为0,0,1,2,n x x nh n =+=L ．要求微分方程数值解，首先要建立数值算法，即对初值问题（1）中的方程离散化，建立求解数值解法的递推公式．一类是计算1n y +时只用到前一点的值n y ，称为单步法；另一类是用到1n y +前面k 点的值11,,,n n n k y y y --+L 称为k 步法.对初值问题（1）式的单步法可用一般形式表示为11(,,,)n n n n n y y h x y y h ϕ++=+⋅，其中多元函数ϕ与(),f x y 有关，当ϕ含有1n y +时，方法是隐式的；若ϕ中不含1n y +，则为显式方法，所以显式单步法可表示为1(,,)n n n n y y h x y h ϕ+=+⋅.（2）设()y x 是初值问题（1）的准确解，称()()()11(,,)n n n n n T y x y x h x y x h ϕ++=--⋅为显式单步法（2）的局部截断误差. 若存在最大正整数p ,使显式单步法（2）式的局部截断误差满足()()()()11,,p n T y x h y x h x y h O h ϕ++=+--=，则称（2）式有p 阶精度.1.2几种常用的数值解法及其分析、比较1.2.1欧拉法与后退欧拉法1）欧拉法：欧拉曾简单地用差分代替微分，即利用公式将初值问题（1）离散化，则问题（1）可化为1(,),n n n n y y h f x y +=+⋅0n x x n h =+⋅， （3）此方法称为欧拉法.欧拉方法的几何意义在数值计算公式中体现了出来.在xy 平面上，一阶微分方程的解()y y x =称作它的积分曲线.积分曲线上一点(),x y 的切线斜率等于函数(),f x y ，按函数(),f x y 在xy 平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.基于上述几何解释，从初始点000(,)P x y 出发，先依方向场在该点的方向上推进到1x x =上一点1P ，再从1P 依方向场的方向推进到2x x =上一点2P ，循环前进便作出一条折线012P PP L ，因此欧拉方法又称为折线法.若初值0y 已知，则由（3）式可逐步算出为了分析计算公式的精确度，通常可用泰勒展开将()1n y x +在n x 处展开，则有()()()()()()2''11,,.2++'=+=++∈n n n n n n n n h y x y x h y x y x h y x x ξξ 在()n n y y x =的前提下，()()()(),,.n n n n n f x y f x y x y x '==可得欧拉法（3）的误差为 容易看出，欧拉法（3）式具有一阶精度.2）向后欧拉方法：如果对微分方程（1）从n x 到1n x +积分，得()()()()11,n n x n n x y x y x f t y t dx ++=+⎰，（4） 如果（4）式右端积分用右矩形公式()()11,n n h f x y x ++⋅近似，则得到另一个公式 ()111,n n n n y y hf x y +++=+， （5） 称为后退欧拉法.值得一提的是:后退欧拉法与欧拉公式有着本质的区别，后者是关于1n y +的直接计算公式，它是显式的，而（5）式的右端含有关于1n y +的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化.具体迭代过程如下：首先利用欧拉公式(0)1(,)+=+⋅n n n n y y h f x y 给出迭代初值(0)1+n y ，把它代入（5）式的右端，使之转化为显式，直接计算得11(1)(0)1(,)+++=+⋅n n n n y y h f x y .如此反复进行，得 (1)()111(,)++++=+⋅k k n n n n y y h f x y 0,1,k =L ，则得到后退欧拉法的迭代公式(0)1(1)()111(,)(,)+++++⎧=+⋅⎨=+⋅⎩n n n n k k n n n n y y h f x y y y h f x y ， 可以看出，后退欧拉法具有一阶精度，且计算比较麻烦.1.2.2梯形方法为得到比欧拉法精确度高的计算公式，在等式（4）式右端积分中若用梯形求积公式近似，并用n y 代替()n y x ，1n y +代替()1n y x +，则得()()111,,2n n n n n n h y y f x y f x y +++=++⎡⎤⎣⎦， （6） 称其为梯形方法.梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为 ()()(0)1(1)()111(,),,2+++++⎧=+⋅⎪⎨⎡⎤=++⎪⎣⎦⎩n n n n k k n n n n n n y y h f x y h y y f x y f x y . （7） 为了分析梯形公式的收敛性，将（6）与（7）式相减，得()()(1)()111111,,2k k n n n n n n h y y f x y f x y +++++++⎡⎤-=-⎣⎦，0,1,2,k =L 因为(),f x y 满足Lipschitz 条件，于是有(1)()11112+++++-≤-k k n n n n hL y y y y ，其中L 为(),f x y 关于y 的Lipschitz 常数.如果选取h 充分小，使得12hL <，则当k →∞时有(1)11+++→k n n y y ，这说明迭代过程（7）式是收敛的[4].容易推导得出梯形法（7）式是二阶方法.经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的.从上述的迭代公式可以看出，每迭代一次都要重新计算(),f x y 的值，而且迭代又要进行若干次，计算相当的复杂.为此，有没有比较简便的计算方法呢？下面给出改进的欧拉方法.1.2.3改进的欧拉方法由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方法是只计算一两次就转入下一步的计算，先用欧拉公式（3）求得一个初步的近似解1n y +，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常称为改进的欧拉公式.具体公式如下()()()11,,,2n n n n n n n n h y y f x y f x y hf x y ++⎡⎤=+++⎣⎦ （8）改进的欧拉法与梯形法一样，是二阶方法.1.2.4 Runge-Kutta 方法由前面讨论可知，从（4）式可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积积点，为此将（4）式的右端用求积公式表示为()()()()11,,n n rx i n i n i x i f x y x dx h c f x h y x h λλ+=≈++∑⎰， （9） 一般来说，点数r 越多，精度越高，上式右端相当于增量函数(),,x y h ϕ，为得到便于计算的显式方法，将公式（9）表示为：()1,,,n n n n y y h x y h ϕ+=+（10） 其中()()1111,,,,,2,r n n i i i n n i i n i n ij j j x y h c K K f x y K f x h y h K i r ϕλμ=-=⎧⎪=⎪⎪=⎨⎪⎛⎫⎪=++= ⎪⎪⎝⎭⎩∑∑L （11） 这里,,i i ij c λμ均为常数. i c 为加权因子，i K 为第i 段斜率，共有r 段.我们把（10）和（11）称为r 级显式Runge-Kutta 法，简称为R-K 方法.下面给出其中最经典最常用的一个公式：()()()11234121324322,6,,,22,,22,.n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎛⎫⎪=++ ⎪⎨⎝⎭⎪⎪⎛⎫=++ ⎪⎪⎝⎭⎪⎪=++⎩（12） Runge-Kutta 方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算1n y +，只用到前面一步的值n y 即可，因此每步的步长可以独立取定.常用的Runge-Kutta 方法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，步长h 可取得大些，求解区间上的总步数可以少些.但Runge-Kutta 方法也有些缺点，比如四阶Runge-Kutta 方法每算一步需要四次计算(),f x y 的值，计算量较大（对于复杂的(),f x y 而言）． 2 数值方法的应用实例[5-9]例1 对于初值问题()1001y y y '=-⎧⎪⎨=⎪⎩，分别用欧拉法、改进的欧拉法，梯形法求()1y 的近似值.解：易得该方程的解析解()10x y x e -=，()1 4.5400e-005y =，为比较，将按不同数值计算方法所得结果列表如下：表 1 三种不同方法的数值结果图 1 三种不同方法数值解与精确解的误差曲线从表1中可以看出：当0.2h=时，三种方法均不稳定，计算结果严重偏离精确值；h=时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；当0.01h≤时，三种0.1方法均稳定，但精确度有区别．可以看出，h越小，计算结果越好，要想计算结果充分接近于解析解还须取较小的h值．图1反映的步长0.01h=时，三种数值方法的所得数值解与解析解在[0,1]区间的误差曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；改进的欧拉法和梯形法精确度都明显高于欧拉法．例2用欧拉法、改进的欧拉法和Runge-Kutta法求解初值问题并比较三种方法的结果．解：方程为1n=-的伯努利方程，可求得解析解为现用MATLAB软件编程，用题目要求的方法求解，可得如下图示结果：图2 （a）步长为0.2时R-K法和解析解比较图2 （b）步长为0.2时改进的Euler法和解析解比较图2 （c）步长为0.2时欧拉法和解析解比较上图2(a)，(b)，(c)描述的是步长为0.2时，用欧拉法、改进的欧拉法，Runge-Kutta 法求解方程所得的数值解与解析解之间的对比图．由图可知，Runge-Kutta法所得数值解曲线和解析解曲线吻合的很好，改进的欧拉法和欧拉法随着计算的进行，数值解和解析解之间误差逐步增大，但改进的欧拉法效果要好于欧拉法．图3 (a) 步长为0.1时Euler法和解析解比较图3 (b) 步长为0.1时改进的Euler法和解析解比较图3 (c) 步长为0.1时Runge-Kutta法和解析解比较上图3 (a)，(b)，(c)描述的是步长为0.1时，用欧拉法、改进的欧拉法，Runge-Kutta法求解方程所得的数值解与解析解之间的对比图．由图可知，改进的欧拉法和Runge-Kutta法所得数值解曲线和解析解曲线吻合的很好，而欧拉法随着计算的进行数值解和解析解之间误差逐步增大．相应的程序如下：主程序x=0:0.2:1;jxj=exp(2*x).*(1./exp(4*x) + (2*x)./exp(4*x)).^(1/2);y=Euler(@ff,0,1,0.2,1);gy=geuler(@ff,0,1,0.2,1);Ry=RK(@ff,0,1,0.2,1); figure(1);plot(x,jxj,x,Ry,'*');figure(2);plot(x,jxj,x,gy,'*');figure(3);p lot(x,jxj,x,y,'*')欧拉法程序function y=Euler(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:ny(i+1)=y(i)+h*feval(f,x(i),y(i));end改进的欧拉法程序function gy=geuler(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:nyp=y(i)+h*feval(f,x(i),y(i));yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=(yp+yc)/2;endgy=y;Runge-Kutta法程序function Ry=RK(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:nk1=feval(f,x(i),y(i));k2=feval(f,x(i)+h/2,y(i)+h*k1/2);k3=feval(f,x(i)+h/2,y(i)+h*k2/2);k4=feval(f,x(i+1),y(i)+h*k3);y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;endRy=y;3 微分方程数值解法在实际生活中的应用3.1应用实例:耐用消费新产品的销售规律模型一种新产品进入市场以后，常常会经历销售量首先慢慢增加然后逐渐慢慢下降的一个过程，由此给出的时间—销售坐标系下的曲线称为产品的生命曲线,它的形状是钟形的.不过对于较耐用的消费品，情况会有所不一样，它的生命曲线会在开始有个小小的高峰，之后是段比较平坦的曲线，先下降，之后再上升，然后达到高峰，因此它是双峰型的曲线.如何理解这种和传统产品的生命曲线理论相冲突的现象？澳大利亚学者斯蒂芬斯与莫赛观察到的购买耐用消费产品的人大概可分为两类：其一是非常善于接受新的事物，称其为“创新型”消费者，他们会经常从产品广告，制造商给出产品的说明书与商店样品中了解了产品功能与性能之后，再决定其否购买；其二是消费者相对保守，称其为“模仿型”消费者，他们往往会根据大部分已经购买产品的消费者实际使用的经验而提供的信息来决定是否购买产品，下面的例子经过分析，建立相应的数学模型，对这种现象给出了科学解释.3.1.1 模型的建立将顾客获得信息大致可分成两类，其一称之为“搜集型”，源自于产告说明、广告，“创新型”顾客获得这类消息后就可作出其是否购买；第二类称之为“体验型”，即消费者使用之后会获得真实体验，常常以口头的形式散播，“模仿型”类顾客会在获得这种信息之后才能决定是否购买.设K 是潜在用户的总数，1K 与2K 分别为 “创新型”与“模仿型”的人数.再设()N t 是时刻t 已经购买的商品顾客的人数，而()1N t 与()2N t 分别为“创新型”与“模仿型”的顾客的人数,再设()1A t 是时刻t 中已获得“搜集型”信息的人数，由于该部分的信息可直接由外部获得，或者可从已获得该类信息的人群中获得，因此类似巴斯模型，从而建立如下方程：()()()()()()11112112,00,,0dA t K A t A t A dtαααα=-+=> ， 获得“搜集型”信息的“创新型”消费者决定其是否购买的行为，满足如下方程： 对于“模仿型”的顾客，可从已经购买该产品“创新型”或者“模仿型”的顾客中获取信息，于是有在这里，忽略消费者购买产品后需一段短暂的使用后才会散播体验信息的滞后作用． 综上，斯蒂芬斯—莫赛模型是一常微分方程组的初值问题：而()()()12N t N t N t =+为时刻t 购买该商品的总人数.3.1.2 模型的求解对于斯蒂芬斯—莫赛模型中()2N t 的解析解则不能求出，于是可以用四阶Runge Kutta -公式求得，且在它的精度要求达到很高情形下求出()2N t .用MATLAB 软件求解，相应的程序及结果如下function RK=RKFC(fc,a,b,h,y0)n=(b-a)/h;x=a:h:b;m=length(y0);y=zeros(n+1,m);y(1,:)=y0;for i=1:nk1=feval(fc,x(i),y(i,:));k2=feval(fc,x(i)+h/2,y(i,:)+h*k1/2);k3=feval(fc,x(i)+h/2,y(i,:)+h*k2/2);k4=feval(fc,x(i+1),y(i,:)+h*k3); y(i+1,:)=y(i,:)+h*(k1+2*k2+2*k3+k4)/6;endRK=y;function f=FC(x,y)k1=50; k2=70;al=0.0013;be=0.0013;ga=0.0015;f=[(k1-y(1))*(al+be*y(1)),ga*(k2-y(2))*(y(1)+y(2))];x=0:0.3:24;RK=RKFC(@FC,0,24,0.3,[0,0])figure(1)；plot(x,RK(:,1),'+',x,RK(:,2),'*')；legend('N1(t)','N2(t)',2) figure(2)；plot(x,RK(:,1)+RK(:,2),'+')；legend('N(t)',2)图 4 ()()12,N t N t 与时间关系图图 5 ()N t 与[0,25]时间段关系图 由此例可以看出，微分方程数值解在实际生活有着广泛的应用，而数值解法中的Runge-Kutta 方法是一种很重要且应用很广泛的算法.结语微分方程数值解是求解微分方程的一种很重要且应用范围很广的方法，而常用的数值解法如欧拉法、改进的欧拉法、梯形法和Runge-Kutta 方法在一定程度上都有自己的优缺点，理解和掌握各种方法的应用范围，用MATLAB 编制各种方法求解实际问题的通用程序，对用微分方程数值解理论解决现实生活中的实际问题有很重要的意义．参考文献[1] 李庆扬,王能超,易大义.数值分析[M].北京:清华大学出版社,清华大学出版社,2001.[2] 冯康.数值计算方法[M].杭州:浙江大学出版社,2003.[3] 封建湖.计算方法典型题分析解集[M].西安:西北工业大学出版社,2000.[4] 胡建伟,汤怀民.微分方程数值解法(第二版)[M].北京:科学出版社,2007.[5] 王能超.计算方法简明教程[M].北京:高等教育出版社,2004.[6] Nash S G.A history of scientific computing.[M] New York:ACM Press,1990.[7] 于丽妮. ODE 问题解析解及数值解的matlab 实现[J]. 电脑知识与技术.2012,8(14):3303-3305.[8] 霍晓成. 常微分方程数值解法的研究[J]. 临沂师范学院学报. 2011,(6):19-23.[6] 王国立, 陈 瑛.非线性微分方程迭代算法及其在物理学中的应用[J]. 长春师范学院学报( 自然科学版). 2006, 25(2):10-12.。
