0计算方法及MATLAB实现简明讲义课件PPS8-1欧拉龙格法
龙哥库塔法or欧拉法求解微分方程matlab实现
龙哥库塔法or欧拉法求解微分⽅程matlab实现举例：分别⽤欧拉法和龙哥库塔法求解下⾯的微分⽅程我们知道的欧拉法（Euler）"思想是⽤先前的差商近似代替倒数"，直⽩⼀些的编程说法即：f(i+1)=f(i)+h*f(x,y)其中h是设定的迭代步长，若精度要求不⾼，⼀般可取0.01。
在定义区间内迭代求解即可。
龙哥库塔法⼀般⽤于⾼精度的求解，即⾼阶精度的改进欧拉法，常⽤的是四阶龙哥库塔，编程语⾔如下：y(i+1)=y(i)+h*(k1+2*K2+2*k3+k4)/6;k1=f(xi,yi)k2=f(xi+h/2,yi+h*k1/2);k3=f(xi+h/2,yi+h*k2/2);k4=f(xi+h,yi+h*k3);设置终⽌条件迭代求解。
matlab实现程序如下：%% 龙哥库塔or欧拉法求解微分⽅程t=0:0.01:3; %⾃变量范围f = [];g = [];f(1) = 0.1; %f初值g(1) = 1; %g初值h=0.001;for i=1:length(t)%% 欧拉法% f(i+1) =f(i)+h*(exp(f(i)*sin(t(i)))+g(i));% g(i+1) =g(i)+h*(exp(g(i)*cos(t(i)))+f(i));%% 龙哥库塔kf1 = exp(f(i)*sin(t(i)))+g(i);%g(i)相当于常值kf2 = exp((f(i)+kf1*h/2)*sin(t(i)+h/2))+g(i);kf3 = exp((f(i)+kf2*h/2)*sin(t(i)+h/2))+g(i);kf4 = exp((f(i)+kf3*h)*sin(t(i)+h))+g(i);f(i+1) = f(i) + h*(kf1+2*kf2+2*kf3+kf4)/6;kg1 = exp(f(i)*cos(t(i)))+f(i);%f(i)相当于常值kg2 = exp((f(i)+kg1*h/2)*cos(t(i)+h/2))+f(i);kg3 = exp((f(i)+kg2*h/2)*cos(t(i)+h/2))+f(i);kg4 = exp((f(i)+kg3*h)*cos(t(i)+h))+f(i);g(i+1) = g(i) + h*(kg1+2*kg2+2*kg3+kg4)/6;endfigure(1)plot(t,f(1:length(t)),t,g(1:length(t)));legend('f函数','g函数')。
matlab教程ppt(完整版)
展示部分与整体的关系，通过扇形面积或角度表 示占比。
三维图形
01
02
03
04
三维散点图
在三维空间中展示两个变量之 间的关系，通过点的位置展示
数据。
三维曲面图
通过曲面表示两个或多个变量 之间的关系，可以展示数据的
分布和趋势。
三维等高线图
表示三维空间中数据的分布和 变化，通过等高线的形状和密
集程度展示数据。
处理运行过程中出现的错误和 异常情况。
通过优化算法和代码结构，提 高程序的运行效率。
对代码进行重新组织，使其更 易于阅读和维护。
03
MATLAB可视化
绘图基础
散点图
描述两个变量之间的关系，通过点的分布展示数 据。
条形图
比较不同类别的数据大小，通过条形的长度或高 度进行比较。
折线图
展示时间序列数据或多个变量之间的关系，通过 线条的走势呈现数据变化。
控制系统仿真
使用MATLAB进行控制系统仿真 ，模拟系统动态性能。
控制系统优化
对控制系统进行优化设计，如权 重优化、多目标优化等。
THANK YOU
感谢聆听
对图像进行几何变换，如缩放、旋转、平移 等操作。
动画制作
帧动画
通过一系列静态图像的连续播放，形 成动态效果。
路径动画
让对象沿指定路径移动，形成动态效 果。
变形动画
让对象从一个形状逐渐变形为另一个 形状，形成动态效果。
交互式动画
允许用户通过交互操作控制动画的播 放、暂停、回放等操作。
04
MATLAB在科学计算中的应用
对函数进行数值积分和微分， 用于解决定积分和微分方程问 题。
数值优化
matlab教程ppt(完整版)
数据处理
应用MATLAB的信号处理和统计 分析函数库，进行数据预处理、
特征提取和模型训练。
机器学习与深度学习
机器学习
介绍MATLAB中的各种机器学习算法，如线性回归、决策 树、支持向量机等，以及如何应用它们进行分类、回归和 聚类。
深度学习
介绍深度学习框架和网络结构，如卷积神经网络（CNN） 、循环神经网络（RNN）等，以及如何使用MATLBiblioteka B进行 训练和部署。感谢观看
THANKS
符号微积分
进行符号微分和积分运算，如极限、导数和 积分。
符号方程求解
使用solve函数求解符号方程。
符号矩阵运算
进行符号矩阵的乘法、转置等运算。
05
MATLAB应用实例
数据分析与可视化
数据分析
使用MATLAB进行数据导入、清 洗、处理和分析，包括描述性统
计、可视化、假设检验等。
可视化
利用MATLAB的图形和可视化工 具，如散点图、柱状图、3D图等
数值求和与求积
演示如何对数值进行求和与求积 操作。
数值计算函数
介绍常用数值计算函数，如sin、 cos、tan等。
方程求解
演示如何求解线性方程和非线性方 程。
03
MATLAB编程基础
控制流
01
02
03
04
顺序结构
按照代码的先后顺序执行，是 最基本的程序结构。
选择结构
通过if语句实现，根据条件判 断执行不同的代码块。
数据分析
数值计算
MATLAB提供了强大的数据分析工具，支 持多种统计分析方法，可以帮助用户进行 数据挖掘和预测分析。
MATLAB可以进行高效的数值计算，支持 多种数值计算方法，包括线性代数、微积 分、微分方程等。
matlab教程ppt完整版
进行图像的裁剪、缩放、旋转等基本操作，以满 足图像处理的需求。
图像处理特效
应用滤波、边缘检测、色彩空间转换等图像处理 技术，提升图像质量或提取图像特征。
程序设计与优化
05
M文件编程基础
M文件概述
01
M文件是MATLAB中用于存储代码和数据的文本文件，具有.m
扩展名。
脚本文件与函数文件
稀疏矩阵压缩
通过压缩存储方式节省内存空间。
稀疏矩阵运算
支持基本的四则运算和矩阵函数。
稀疏矩阵应用
在数值计算、图像处理等领域有广泛应用。
数值计算与函数分
03
析
多项式运算及函数拟合
多项式表示与运算
介绍如何在MATLAB中创建多项 式、进行多项式四则运算以及多
项式求值。
函数拟合方法
详细阐述最小二乘法、梯度下降法 等函数拟合方法，并给出相应的 MATLAB实现代码。
使用plot3、mesh、surf等函数 绘制三维曲线、曲面图。
三维图形视角调整
通过view、rotate等函数调整三 维图形的观察角度，以便更好地
展示数据特征。
三维图形样式设置
设置颜色映射、透明度、光照效 果等，提升三维图形的视觉效果
。
特殊图形绘制技巧
极坐标与对数坐标绘图
使用polar、semilogx、semilogy等函数绘制极坐标图和对数坐 标图，适应不同类型的数据展示需求。
使用`dsolve`命令求解常微分方程，使用 `pdepe`等命令求解偏微分方程，分析物理 现象和工程问题。
MATLAB高级功能
07
与应用
MATLAB编译器使用指南
MATLAB编译器介绍
matlab教程ppt(完整版)
矩阵减法：两个相同大小 的矩阵可以进行减法运算 ，例如D=A-B。
矩阵的分解与特征值
详细描述
矩阵分解：将一个复杂的矩阵分 解为几个简单的、易于处理的矩 阵，例如LU分解、QR分解等。
特征值：矩阵的特征值是该矩阵 的一个重要的数值属性，可以用 于分析矩阵的性质和特征。
矩阵运算
介绍矩阵的创建、索引、算术 运算和逻辑运算等操作。
控制流
介绍if语句、for循环和while 循环等控制流结构的使用方法 。
02
MATLAB编程
变量与数据类型
01
02
03
变量命名规则
MATLAB中的变量名以字 母开头，可以包含字母、 数字和下划线，但不能包 含空格。
数据类型
MATLAB支持多种数据类 型，如数值型、字符型、 逻辑型和单元数组等。
matlab教程PPT(完整版)
汇报人：可编辑 2023-12-26
目 录
• MATLAB基础 • MATLAB编程 • MATLAB矩阵运算 • MATLAB图像处理 • MATLAB数值分析 • MATLAB应用实例
01
MATLAB基础
MATLAB简介
MATLAB定义
MATLAB应用领域
MATLAB是一种用于算法开发、数据 可视化、数据分析和数值计算的编程 语言和环境。
函数编写
01
02
03
04
函数定义
使用`function`关键字定义函 数，指定输入输出参数。
函数体
在函数定义中编写实现特定功 能的代码。
函数调用
通过函数名和输入参数调用自 定义函数。
欧拉方法matlab
欧拉方法matlab
欧拉方法是一种数值解微分方程的方法，适用于一阶常微分方程和某些高阶微分方程。
本文将介绍如何使用MATLAB实现欧拉方法。
首先，我们需要定义微分方程和初始条件。
例如，对于一阶常微分方程dy/dx = f(x,y)，初始条件为y(x0) = y0，我们可以定义函数f和初始值x0和y0：
function dydx = f(x,y)
dydx = ... % 定义f(x,y)的具体表达式
end
x0 = ... % 初始值x0
y0 = ... % 初始值y0
接下来，我们需要设置步长和计算区间。
步长越小，计算结果越精确，但计算量也会增加。
计算区间可以根据需要设置。
h = ... % 步长
xspan = ... % 计算区间
然后，我们可以使用欧拉方法计算微分方程的数值解。
欧拉方法的基本思想是在每个步长上使用当前值和导数的近似值（如斜率）来估计下一个值。
首先，我们可以定义一个数组来存储计算结果，设初始值为y0： y = [y0];
然后，我们可以使用for循环来计算每个步长上的结果。
在每个步长上，我们需要计算当前值的导数，以及使用欧拉方法计算下一个
值。
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);```。
matlab实例讲解欧拉法求解微分方程
欧拉法是数值分析中常用的一种方法，用于求解常微分方程的数值解。
在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中使用欧拉法求解微分方程是非常简单而直观的。
解微分方程欧拉法RK法及其MATLAB实例
解微分方程的欧拉法，龙格-库塔法及其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教程ppt(完整版) (3)
数值积分与微分
数值积分
使用MATLAB的`integral`函数进 行数值积分，可以选择不同的积
分方法。
数值微分
可以使用差分法或`diff`函数进行 数值微分。
符号积分与微分
使用符号计算工具箱中的函数， 如`syms`、`int`和`diff`，进行符
号积分和微分。
常微分方程求解
欧拉法
简单的一阶常微分方程的初值问题可以使用欧拉法求解。
图形可视化
MATLAB具有强大的图形可视化功能，支 持多种图形类型和交互操作。
编程语言
MATLAB是一种高级编程语言，具有丰富 的函数库和工具箱。
数据分析
MATLAB提供了多种数据分析工具，包括 数据导入、处理、分析和可视化。
MATLAB的应用领域
科学计算
广泛应用于数学、物理、工程等 领域。
控制系统设计
短时傅里叶变换
通过在时间上滑动窗口并对每个窗口内的信号进 行傅里叶变换，实现信号的时频分析。
小波变换
利用小波基函数的特性，对信号进行多尺度分析 ，从而在时频域上展示信号的细节。
信号滤波与变换
数字滤波器设计
使用MATLAB中的滤波器设计工具，如butterworth、 chebyshev等，设计数字滤波器以实现信号的滤波。
03 多目标优化
使用`gamultiobj`函数求解多目 标最优化问题。
0 最小二乘问题 4使用`lsqlin`或`lsqnonlin`函数
求解线性或非线性最小二乘问 题。
05
MATLAB在信号处理中的应用
信号的时频分析
信号的时频表示
将信号从时间域转换到时频域，以便更好地理解 和分析信号的特性。
matlab教程ppt(完整版)
• 2002年7月，推出了Matlab 6.5(R13)，在这一版本中Simulink升级到了5.0，性能有 了很大提高，另一大特点是推出了JIT程序加速器，Matlab的计算速度有了明显的 提高。 • 2005年9月，推出了MAILAB 7.1(Release14 SP3)，在这一版本中Simulink升级到了 6.3，软件性能有了新的提高，用户界面更加友好。值得说明的是，Matlab V7.1版 采用了更先进的数学程序库，即“LAPACK”和“BLAS”。
MATLRAeBal-TToiPmorleobcoWexsoessriknsghBolpo是ck一很set种可等实能，时已详代有见码人M生A将T成你LA工要B具做在，的线它应帮能用助够程文根序据作成工具箱了。 MATLS成Aim实BuCl时ino档k应m模p。用i型le程r生序成。程序源代码，并打包、编译所生成的源代码生 Simulink Stateflow从是现基有于的有Si限mu状lin态k 机和理Sta论te针flo对w自复动杂生成C语言程序代码的功能、
解微分方程欧拉法RK法及其MATLAB实例
解微分方程的欧拉法，龙格-库塔法及其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决定。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第8章常微分方程初值问题数值解法8.1 引言8.2 欧拉方法8.3 龙格-库塔方法8.4 单步法的收敛性与稳定性8.5 线性多步法8.1 引 言考虑一阶常微分方程的初值问题00(,),[,],().y f x y x a b y x y '=∈=（1.1） （1.2）如果存在实数 ，使得121212(,)(,).,R f x y f x y L y y y y -≤-∀∈（1.3）则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点)(x y<<<<<+121n n x x x x 上的近似值 .,,,,,121+n n y y y y 相邻两个节点的间距 称为步长.n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ),2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”.即求解过程顺着节点排列的次序一步一步地向前推进.00(,),[,],().y f x y x a b y x y '=∈=描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式.1+n y 一类是计算时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值， 1+n y k 11,,,+--k n n n y y y 称为 步法.k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题.n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式.),(y x f y ='8.2 简单的数值方法8.2.1 欧拉法与后退欧拉法积分曲线上一点 的切线斜率等于函数 的 值.),(y x ),(y x f 如果按函数 在平面上建立一个方向场，那么， ),(y x f xy 积分曲线上每一点的切线方向均与方向场在该点的方向相一致.在 平面上，微分方程的解 称 作它的积分曲线.xy )(x y y =),(y x f y ='基于上述几何解释，从初始点 出发， ),(000y x P 先依切线 在该点的方向推进到 上一点 ，然后再 1x x =1P 从 依切线 的方向推进到 上一点 ，循此前进 1P 2x x =2P 做出一条折线 （图8-1).210P P P 图8-179一般地，设已做出该折线的顶点 ，过 依 切线 的方向再推进到 ，显然两个顶点的坐标有关系 ),(n n n y x P nP ),(111+++n n n y x P 1,+n n P P 11n n n n y y x x ++-=-反解，得).,(1n n n n y x hf y y +=+（2.1）这就是著名的欧拉（Euler ）方法. 欧拉法实际上是对常微分方程中的导数用差商近似，即1()()n n y x y x h+-≈直接得到的.(,)d d n n x y y x =(,)n n f x y 00(,),[,],().y f x y x a b y x y '=∈=1n n y y h +-=()(,())n n n y x f x y x '=111(,)n n n n y y hf x y +=+例欧拉方法13若初值已知，则依公式（2.1）可逐步算出 0y ),,(0001y x hf y y +=),,(1112y x hf y y +=).,(1n n n n y x hf y y +=+（2.1）例1 求解初值问题⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y x y y （2.2）解 欧拉公式的具体形式为).2(1nnn n n y x y h y y -+=+取步长 ，计算结果见表8-1.1.0=h 初值问题（2.2）的解为 ，按这个解析式 子算出的准确值 同近似值 一起列在表8-1中，两者相比较可以看出欧拉方法的精度不高.x y 21+=)(n x y n y 1(,)n n n n y y hf x y +=+12y x=+81()()0.1 1.1000 1.09540.6 1.5090 1.48320.2 1.1918 1.18320.7 1.5803 1.54920.3 1.2774 1.26490.8 1.6498 1.61250.4 1.3582 1.34160.9 1.7178 1.67330.51.43511.41421.01.78481.7321n nn n nn x y y x x y y x -表计算结果对比 还可以通过几何直观来考察欧拉方法的精度.假设 ，即顶点 落在积分曲线 上，那么，按欧拉方法做出的折线 便是 过点 的切线（图8-2）. )(n n x y y =nP )(x y y =1+n n P P )(x y y =nP图8-2从图形上看，这样定出的顶点 显著地偏离了原来 的积分曲线，可见欧拉方法是相当粗糙的.1+n P 为了分析计算公式的精度，通常可用泰勒展开将)(1+n x y在处展开，则有n x)()(1h x y x y n n +=+在 的前提下， )(n n x y y =)())(,(),(n n n n n x y x y x f y x f '==11()n n y x y ++-=称为此方法的局部截断误差.()()n n y x y x h '=++于是可得欧拉法（2.1）的误差 （2.3）),(22n x y h ''≈如果对方程 从 到积分，得 n x 1+n x ),(y x f y =').,(1n n n n y x hf y y +=+（2.1） 2()2n h y ξ''1(,)n n n x x ξ+∈1(,)n n n n y y hf x y +=+(上一值精确相等时，估计下一值误差!) 2()2n h y ξ''1()n y x +=（2.4）1()(,())n nx n xy x f t y t dt++⎰18右端积分用左矩形公式 近似.))(,(n n x y x f h 再以 代替 n y ),(n x y 如果在（2.4）中右端积分用右矩形公式 ))(,(11++n n x y x f h 111(,)n n n n y y hf x y +++=+（2.5）称为后退的欧拉法.代替也得到（2.1）， )(1+n x y 1+n y 局部截断误差也是（2.3）. 近似，则得另一个公式它也可以通过利用均差近似导数 ，即11()()n n n ny x y x x x ++-≈-)(1+'n x y 1(,)n n n n y y hf x y +=+2()2n hy ξ''1()(,())n nx n x y x f t y t dt++⎰111()(,())n n n y x f x y x +++'=欧拉公式是关于 的一个直接的计算公式，这类公式称作是显式的；1+n y 后退欧拉公式的右端含有未知的 ，它是关于 的一个函数方程，这类公式称作是隐式的.1+n y 1+n y 1(,)n n n n y y hf x y +=+111(,)n n n n y y hf x y +++=+隐式方程通常用迭代法求解，而迭代过程的实质是逐步显示化. 设用欧拉公式),()0(1n n n n y x f h y y+=+给出迭代初值 ，用它代入（2.5）式的右端，使之转化 为显式，直接计算得)0(1+n y ),,()0(11)1(1++++=n n n n yx f h y y然后再用 代入（2.5）式，又有)1(1+n y ).,()1(11)2(1++++=n n n n yx f h y y111(,)n n n n y y hf x y +++=+（2.5）21如此反复进行，得).,1,0(),,()(11)1(1=+=++++k yx f h y yk n n n k n （2.6）由于 对 满足李普希茨条件(1.3). ),(y x f y 由(2.6)减(2.5)得(1)11k n n yy +++-=.1)(1++-≤n k n y y hL 由此可知，只要迭代法(2.6)就收敛到解 . 1<hL 1+n y 从积分公式可以看到后退欧拉方法的公式误差与欧拉法是相似的.111(,)n n n n y y hf x y +++=+（2.5） ()1111(,)(,)k n n n n h f x y f x y ++++-121212(,)(,).,Rf x y f x y L y y y y -≤-∀∈8.2.2 梯形方法 若用梯形求积公式近似等式（2.4）右端的积分并分别用 代替 则可得到比欧拉法 精度高的计算公式1,+n n y y ),(),(1+n n x y x y 111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） 称为梯形方法. 梯形方法是隐式单步法，可用迭代法求解.⎰+1))(,(n nx x dtt y t f .))(,()()(11⎰++=+n nx xn n dt t y t f x y x y （2.4） 1()(,())n nx n xy x f t y t dt++⎰111(,())[(,)(,)]2n n x n n n n x hf t y t dt f x y f x y +++≈+⎰23⎪⎪⎩⎪⎪⎨⎧+=+);,()0(1n n n ny x f h y y 为了分析迭代过程的收敛性，将（2.7）与（2.8）式相减，得)],(),([2)(11)1(1k nn n n n k n y x f y x f h y y ++++++=（2.8） ).,2,1,0( =k 同后退的欧拉方法一样，仍用欧拉方法提供迭代初值， 则梯形法的迭代公式为111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） (1)()111111[(,)(,)]2k k n n n n n n h y yf x y f x y +++++++-=-如果选取 充分小，使得h 12hL<则当 时有 ， ∞→k 1)(1++→n k n y y 此时迭代过程是收敛的.于是有(1)11k n n y y+++-≤式中 为关于 的李普希茨常数. ),(y x f y L (1)()111111[(,)(,)]2k k n n n n n n h y yf x y f x y +++++++-=-121212(,)(,).,Rf x y f x y L y y y y -≤-∀∈()112k n n hL y y ++-计算方法及MATLAB 实现8.2.3 改进欧拉公式梯形方法虽然提高了精度，但其算法复杂. 在应用迭代公式（2.8）进行实际计算时，每迭代一次，都要重新计算函数 的值.),(y x f 为了控制计算量，通常只迭代一两次就转入下一步的 计算，这就简化了算法.具体地，先用欧拉公式求得一个初步的近似值 , 1n y + 而迭代又要反复进行若干次，计算量很大，而且往往难以预测.称之为预测值，)],(),([2)(11)1(1k n n n n n k n y x f y x f h y y ++++++=计算方法及MATLAB 实现这样建立的预测-校正系统通常称为改进的欧拉公式： 预测值 的精度可能很差，再用梯形公式（2.7）将它校正一次，即按（2.8）式迭代一次得 ，这个结果称校正值. 1+n y 1+n y 预测 ),,(1n n n n y x f h y y +=+校正)].,(),([2111+++++=n n n n n n y x f y x f hy y （2.9）也可以表为下列平均化形式),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(1y y y +=111[(,)(,)]2n n n n n n h y y f x y f x y +++=++（2.7） )],(),([2)(11)1(1k n n n n n k n y x f y x f h y y ++++++=(参阅教材392页!)计算方法及MATLAB 实现例2 用改进的欧拉方法求解初值问题（2.2）.解 这里 改进的欧拉公式为⎪⎪⎪⎩⎪⎪⎪⎨⎧-+=),2(nnn n p y x y h y y ⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y x y y （2.2）2(,)(),xf x y y y=-),2(1p n p n c y x y h y y +-+=).(211c p n y y y +=+),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(211c p n y y y +=+计算方法及MATLAB 实现282(,)xf x y y y=-),,(n n n p y x f h y y +=),,(1p n n c y x f h y y ++=).(211c p n y y y +=+取 ，计算结果见表8-2.1.0=h 同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度.()()0.11.0959 1.09540.6 1.4860 1.48320.2 1.1841 1.18320.7 1.5525 1.54920.3 1.2662 1.26490.8 1.6153 1.61250.4 1.3434 1.34160.9 1.6782 1.67330.51.41641.41421.01.73791.7321n n n n n n x y y x x y y x -表8212y x=+9.2.4 单步法的局部截断误差与阶 初值问题(1.1)，(1.2)的单步法可用一般形式表示为),,,,(11h y y x h y y n n n n n +++=ϕ（2.10）其中多元函数 与 有关，ϕ),(y x f 当 含有 时，方法是隐式的，若不含 则为显式方法，ϕ1+n y 1+n y ),,,(1h y x h y y n n n n ϕ+=+（2.11） 称为增量函数， ),,(h y x ϕ).,(),,(y x f h y x =ϕ所以显式单步法可表示为 例如对欧拉法（2.1）有 它的局部截断误差已由(2.3)给出.⎩⎨⎧=='.)(),,(00y x y y x f y （1.1） （1.2） ).,(1n n n n y x hf y y +=+（2.1） )(2)(211n n n y h y x y ξ''=-++（2.3） ),(22n x y h ''≈对一般显式单步法则可如下定义. 定义1 设 是初值问题(1.1)，(1.2)的准确解， 称 )(x y y =)),(,()()(11h x y x h x y x y T n n n n n ϕ--=++（2.12） 为显式单步法（2.11）的局部截断误差. 之所以称为局部的，是假设在 前各步没有误差.1+n T n x 当 时，计算一步，则有 )(n n x y y =11()n n y x y ++-=1()()(,(),)n n n n y x y x h x y x h ϕ+=--⎩⎨⎧=='.)(),,(00y x y y x f y （1.1）（1.2） ),,,(1h y x h y y n n n n ϕ+=+（2.11）1()[(,,)]n n n n y x y h x y h ϕ+-+1.n T +=1(,)n n n n y y hf x y +=+).,(),,(y x f h y x =ϕ在前一步精确的情况下用显式单步公式计算产生的公式误差. 根据定义，欧拉法的局部截断误差 ))(,()()(11n n n n n x y x hf x y x y T --=++即为（2.3）的结果. 这里 称为局部截断误差主项. )(22n x y h '')(2h O T =局部截断误差可理解为用显式单步方法计算一步的误差，即 ()()n n y x h y x =+-22()()2n h y x o h ''=+显然 ),,,(1h y x h y y n n n n ϕ+=+（2.11）)(2)(211n n n y h y x y ξ''=-++（2.3） ),(22n x y h ''≈()n hy x '-(,)y f x y '=1()()n n y x y x h +=+(泰勒公式) (小o 高阶，大O 同阶)计算方法及MATLAB 实现 定义2 设 是初值问题(1.1)，(1.2)的准确解， 若存在最大整数使显式单步法(2.11)的局部截断误差满足 )(x y p ),(),,()()(11++=--+=p n h O h y x h x y h x y T ϕ（2.13） 则称显式单步法具有 阶精度. p 若将（2.13）展开式写成),())(,(211++++=p p n n n h O h x y x T ψ则 称为局部截断误差主项. 1))(,(+p n n hx y x ψ 以上定义对隐式单步法（2.10）也是适用的.⎩⎨⎧=='.)(),,(00y x y y x f y （1.1） （1.2） ),,,,(11h y y x h y y n n n n n +++=ϕ（2.10） 1(,,)n n n n y y h x y h ϕ+=+(小o 高阶，大O 同阶)对后退欧拉法（2.5）其局部截断误差为))(,()()(1111++++--=n n n n n x y x hf x y x y T 这里 ，是1阶方法，局部截断误差主项为 . 1=p )(22n x y h ''-23()()()2n n h hy x y x O h '''=++).()(232h O x y h n +''-=),,(111++++=n n n n y x hf y y （2.5）2[()()()]n n h y x hy x O h '''-++对梯形法（2.7）有)]()([2)()(111+++'+'--=n n n n n x y x y h x y x y T 所以梯形方法是二阶的，其局部误差主项为 )(123n x y h '''-23()()()23!n n n h h hy x y x y x ''''''=++).()(1243h O x y h n +'''-=)],,(),([2111+++++=n n n n n n y x f y x f h y y （2.7） 24[()()()()]()22n n n n h h y x y x hy x y x O h '''''''-++++3364h h -8.3龙格-库塔方法8.3.1 显式龙格-库塔法的一般形式 上节给出了显式单步法的表达式),,,(1h y x h y y n n n n ϕ+=+其局部截断误差为),(),,()()(11++=--+=p n h O h y x h x y h x y T ϕ对欧拉法 ，即方法为阶， )(21h O T n =+1=p ))].,(,(),([21n n n n n n n n y x hf y h x f y x f h y y ++++=+（3.1）若用改进欧拉法，它可表为此时增量函数 ))].,(,(),([21),,(n n n n n n n n y x hf y h x f y x f h y x +++=ϕ（3.2）与欧拉法的 相比，增加了计算一个 右函数 的值，可望 .),(),,(n n n n y x f h y x =ϕf 2=p 若要使得到的公式阶数 更大， 就必须包含更多的值.p ϕf ,))(,()()(11⎰+=-+n n x x n n dx x y x f x y x y （3.3）从方程 等价的积分形式（2.4），即),(y x f y ='若要使公式阶数提高，就必须使右端积分的数值求积公式 精度提高，必然要增加求积节点.为此可将（3.3）的右端用求积公式表示为线性组合1(,())n n x x f x y x dx +≈⎰点数 越多，精度越高， r 上式右端相当于增量函数 ， ),,(h y x ϕ为得到便于计算的显式方法，可类似于改进欧拉法，将公式表示为),,,(1h y x h y y n n n n ϕ+=+（3.4）其中1(,())r i n i n i i h c f x h y x h λλ=++∑,))(,()()(11⎰+=-+n n x x n n dx x y x f x y x y （3.3）(可变步长因子!) (组合系数!),),,(1∑==r i i i n n K c h y x ϕ),,(1n n y x f K =),(11∑-=++=i j j ij n i n i K h y h x f K μλ（3.5）,,,2r i =这里 均为常数. ij i i c μλ,, (3.4)和(3.5)称为 级显式龙格-库塔（Runge-Kutta ）法， r 简称R-K 方法.当 时，就是欧拉法，),(),,(,1n n n n y x f h y x r ==ϕ此时方法的阶为. 1=p 当 时，改进欧拉法(3.1)，(3.2)也是其中的一种，2=r ),,,(1h y x h y y n n n n ϕ+=+（3.4） ),,(1h y x h y y n n n n ϕ+=+))].,(,(),([21),,(n n n n n n n n y x hf y h x f y x f h y x +++=ϕ).,(1n n n n y x hf y y+=+1(,,)n n n n y y h x y h ϕ+=+(步长调节因子) (组合系数) (斜率组合系数)8.3.2 二阶显式R-K 方法 1(,,)rn n i ii x y h c K ϕ==∑1(,,)n n n n y y h x y h ϕ+=+ 对 的R-K 方法，计算公式如下2=r 11122122211()(,)(,)n n n n n n y y h c K c K K f x y K f x h y hK λμ+=++⎧⎪=⎨⎪=++⎩（3.6）这里 均为待定常数.21221,,,μλc c ),(11∑-=++=i j j ij n i n i K h y h x f K μλ若取 得计算公式,2/1,1,021221====μλc c 12n n y y hK +=+称为中点公式，相当于数值积分的中矩形公式.上式也可表示为1(,(,))22n n n n n n h hy y hf x y f x y +=+++1(,)n n K f x y =（3.10）21(,)22n n h hK f x y K =++)).,(2,2(1n n n n n n y x f hy h x hf y y +++=+11122122211()(,)(,)n n n n n n y y h c K c K K f x y K f x h y hK λμ+=++⎧⎪=⎨⎪=++⎩继续上述过程，经过较复杂的数学演算，可以导出各 种四阶龙格-库塔公式，下列经典公式是其中常用的一个：可以证明其截断误差为 . )(5h O 四阶龙格-库塔方法的每一步需要计算四次函数值 ，f ⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+).,(),2,2(),2,2(),,(),22(6342312143211hK y h x f K K h y h x f K K h y h x f K y x f K K K K K h y y n n n n n n n n n n （3.13） (凸线性组合!) (中点公式!) (欧拉公式!)例 设取步长 ，从 到 用四阶龙格-库塔方法求解 初值问题 ⎪⎩⎪⎨⎧=<<-='.1)0(),10(2y x y xy y 解 这里 ，经典的四阶龙格-库塔公式为 yx y y x f 2),(-=),22(643211K K K K hy y n n ++++=+,2nx y K -=2.0=h 0=x 1=x 112341213243(22),6(,),(,),22(,),22(,).n n n nn n n n n n h y y K K K K K f x y h h K f x y K h hK f x y K K f x h y hK +⎧=++++⎪⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪⎪=++⎪⎩计算方法及MATLAB 实现,222223K h y h x K h y K n n n ++-+=,222112K h y h x K hy K n n n ++-+=.)(2334hK y h x hK y K n n n ++-+= 表8-3列出了计算结果，同时了出了相应的精确解. 龙格-库塔方法的精度较高.112341213243(22)6(,),(,),22(,),22(,).n n n nn 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 +⎧=++++⎪⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪⎪=++⎪⎩yxy y x f 2),(-=,21nnn y x y K -=83()0.2 1.1832 1.18320.4 1.3417 1.34160.6 1.4833 1.48320.8 1.6125 1.61251.01.73211.7321n n n x y y x 表计算结果。