二阶常微分方程的数值求解
二阶常微分方程的数值求解 ppt课件
zk 1
zk
h 6
( L1
2 L2
2 L3
L4 ),
K1 zk , L1 f ( xk , yk , zk ),
K
2
zk
h 2
L1 ,
h
h
h
L2 f ( xk 2 , yk 2 K1 , zk 2 L1 ),
Байду номын сангаас
K3
zk
h 2
L2 ,
h
h
h
L3 f ( xk 2 , yk 2 K 2 , zk 2 L2 ),
%***** R-K4 Method*****% K1=RK_z(i); L1=rightf_sy二s1阶(x常(i微),分R方K程_y的(数i),值R求K解_z(i)); K2=RK_z(i)+0.5*h*L1;
% K1 and 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 L2
clc;clear; h=0.1; a=0;b=2; x=a:h:b; y(1)=1; z(1)=-1; for i=1:length(x)-1
y(i+1)=y(i)+h*z(i); z(i+1)=z(i)+h*y(i); end plot(x,y,'r+',x,exp(-x),'k-'); xlabel('Variable x'); ylabel('Variable y');
二阶常微分方程的数值求解
二阶常微分方程的数值求解
一. 教学要求
求常微分方程的数值解
求常微分方程的数值解一、背景介绍常微分方程（Ordinary Differential Equation，ODE）是描述自然界中变化的数学模型。
常微分方程的解析解往往难以求得，因此需要寻找数值解来近似地描述其行为。
求解常微分方程的数值方法主要有欧拉法、改进欧拉法、龙格-库塔法等。
二、数值方法1. 欧拉法欧拉法是最简单的求解常微分方程的数值方法之一。
它基于导数的定义，将微分方程转化为差分方程，通过迭代计算得到近似解。
欧拉法的公式如下：$$y_{n+1}=y_n+f(t_n,y_n)\Delta t$$其中，$y_n$表示第$n$个时间步长处的函数值，$f(t_n,y_n)$表示在$(t_n,y_n)$处的导数，$\Delta t$表示时间步长。
欧拉法具有易于实现和理解的优点，但精度较低。
2. 改进欧拉法（Heun方法）改进欧拉法又称Heun方法或两步龙格-库塔方法，是对欧拉法进行了精度上提升后得到的一种方法。
它利用两个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。
改进欧拉法的公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\Delta t,y_n+k_1\Delta t)$$$$y_{n+1}=y_n+\frac{1}{2}(k_1+k_2)\Delta t$$改进欧拉法比欧拉法精度更高，但计算量也更大。
3. 龙格-库塔法（RK4方法）龙格-库塔法是求解常微分方程中最常用的数值方法之一。
它通过计算多个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。
RK4方法是龙格-库塔法中最常用的一种方法，其公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_1\Delta t}{2})$$ $$k_3=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_2\Delta t}{2})$$ $$k_4=f(t_n+\Delta t,y_n+k_3\Delta t)$$$$y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\Delta t$$三、数值求解步骤对于给定的常微分方程，可以通过以下步骤求解其数值解：1. 确定初值条件：确定$t=0$时刻的函数值$y(0)$。
微分方程数值解使用数值方法求解微分方程
微分方程数值解使用数值方法求解微分方程微分方程是描述自然现象中变化的数学模型，它是数学和科学研究中的重要工具。
然而，许多微分方程并没有精确的解析解，因此需要使用数值方法来近似求解。
本文将介绍一些常用的数值方法来求解微分方程，包括欧拉方法、改进的欧拉方法和龙格-库塔方法。
一、欧拉方法欧拉方法是最简单、最基础的数值方法之一。
它基于微分方程解的定义，通过离散化自变量和因变量来逼近解析解。
假设我们要求解的微分方程为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```三、龙格-库塔方法龙格-库塔方法是一种更加精确的数值方法，它通过计算多个递推式的加权平均值来逼近解析解。
常微分方程数值解法2线性多步法
03
常见的线性多步法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最简单的一种方法，它基于线性近似，通过已知的函 数值来估计新的函数值。
详细描述
欧拉方法的基本思想是利用已知的函数值来估计下一个点的函数值。具体来说，假设我 们有一个函数 (y = f(x))，在已知 (x_0) 处的函数值 (y_0 = f(x_0)) 的情况下，欧拉方法 通过线性插值来估计 (x_1) 处的函数值 (y_1)，即 (y_1 = y_0 + h f(x_0))，其中 (h) 是
05
线性多步法的优缺点
优点
稳定性好
线性多步法在处理常微分方程时具有较好的数值稳定性， 能够有效地抑制数值振荡，提高计算结果的精度。
01
易于实现
线性多步法的计算过程相对简单，易于 编程实现，适合于大规模数值计算。
02
03
精度可调
通过选择不同的步长和线性多步法公 式，可以灵活地调整计算结果的精度， 满足不同的数值模拟需求。
改进方法的收敛性
研究收敛性条件
深入研究线性多步法的收敛性条件，了解哪些情况下方法可能不收 敛，并寻找改进措施。
优化迭代算法
通过优化迭代算法，提高方法的收敛速度和精度，减少迭代次数， 提高计算效率。
引入预处理技术
利用预处理技术对线性系统进行预处理，改善系统的条件数，提高方 法的收敛性。
拓展应用领域
在工程问题中的应用
控制系统设计
在工程领域中，线性多步法可以用于控制系统设计，通过 建立控制系统的数学模型，设计控制算法和控制器，实现 系统的稳定性和性能优化。
常微分方程数值解
常微分方程数值解常微分方程数值解是数学中的一门重要学科，主要研究如何求解常微分方程，在科学计算中有着重要的应用。
常微分方程模型是自然界中广泛存在的现象描述方法，有着广泛的应用领域。
比如，在物理学中，运动中的物体的位置、速度和加速度随时间的关系就可以通过微分方程描述；在经济学中，经济变化随时间的变化也可以用微分方程来描述。
而常微分方程数值解的求解方法则提供了一种快速、高效的计算手段。
一、常微分方程数值解的基本概念常微分方程就是一个描述自变量（通常是时间）与其导数之间关系的方程。
其一般形式如下：$\frac{dy}{dt} = f(y,t)$其中 $f(y,t)$ 是一个已知的函数。
常微分方程数值解就是对于一个常微分方程，对其进行数字计算求解的方法。
常微分方程数值解常使用数值积分的方法来求解。
由于常微分方程很少有解析解，因此数值解的求解方法显得尤为重要。
二、常微分方程数值解的求解方法常微分方程数值解的求解方法很多，以下介绍其中两种方法。
1.欧拉法欧拉法是最简单的一种数值算法，其思想是通过将一个微分方程转化为一个数值积分方程来求解。
其数值积分方程为：$y_{i+1}=y_i+hf(y_i,t_i)$其中 $h$ 为步长，可以理解为每次计算的间隔。
欧拉法的主要缺点是其精度比较低，收敛速度比较慢。
因此，当需要高精度的数值解时就需要使用其他的算法。
2.级数展开方法级数展开法是通过将一个待求解的微分方程进行Taylor级数展开来求解。
通过对Taylor级数展开的前若干项进行求和，可以得到微分方程与其解的近似解。
由于级数展开法的收敛速度很快，因此可以得到相对较高精度的数值解。
但是，当级数过多时，会出现截断误差。
因此，在实际应用中需要根据所需精度和计算资源的限制来选择适当的级数。
三、常微分方程数值解的应用常微分方程数值解在现代科学技术中有着广泛的应用。
以下介绍其中两个应用领域。
1.物理建模常微分方程的物理建模是常见的应用领域。
二阶代数精度的证明
二阶代数精度的证明一、引言在数值计算中，代数精度是评估数值方法误差的一个重要指标。
对于一个数值方法而言，如果它能够以二阶代数精度逼近解析解，那么我们可以说该方法具有较高的精度。
本文将探讨二阶代数精度的证明，包括定义、推导和具体例子。
二、定义2.1 代数精度代数精度是指数值方法逼近解析解的误差大小。
对于一个数值方法，如果它以二阶代数精度逼近解析解，那么它的误差可以表示为：|E(ℎ)|=C⋅ℎ2其中，E(ℎ)是数值方法的误差，ℎ是步长，C是一个常数。
2.2 步长步长ℎ是指在数值方法中，自变量的取值间隔。
对于微分方程求解问题，步长即是网格点之间的距离。
较小的步长可以提高数值方法的精度，但也会增加计算量。
三、推导3.1 泰勒展开为了证明一个数值方法具有二阶代数精度，我们需要使用泰勒展开。
泰勒展开可以将一个函数在某一点的附近用多项式表示，具体形式如下：f(x+ℎ)=f(x)+ℎf′(x)+ℎ22f″(x)+O(ℎ3)其中，f′(x)和f″(x)分别表示函数f(x)在点x处的一阶和二阶导数。
3.2 数值方法的逼近误差假设我们有一个数值方法，用于求解微分方程。
我们将步长ℎ取得很小，那么在每个网格点处，数值方法的逼近误差可以表示为：E i(ℎ)=u(x i)−u i其中，u(x i)是解析解在网格点x i处的值，u i是数值方法在网格点x i处的值。
3.3 利用泰勒展开推导逼近误差我们可以利用泰勒展开推导数值方法的逼近误差。
以一阶导数为例，将f(x+ℎ)展开到二阶项，可以得到：f(x+ℎ)=f(x)+ℎf′(x)+ℎ22f″(x)+O(ℎ3)将上式中的f(x)替换为u(x)，f′(x)替换为u′(x)，f″(x)替换为u″(x)，得到：u(x+ℎ)=u(x)+ℎu′(x)+ℎ22u″(x)+O(ℎ3)将上式中的x替换为x i，ℎ替换为ℎi，得到：u(x i+ℎi)=u(x i)+ℎi u′(x i)+ℎi22u″(x i)+O(ℎi3)将上式左边的u(x i+ℎi)替换为u i，得到：u i=u(x i)+ℎi u′(x i)+ℎi22u″(x i)+O(ℎi3)将上式两边的u(x i)相减，得到：E i(ℎi)=ℎi u′(x i)+ℎi22u″(x i)+O(ℎi3)将上式整理，得到：E i(ℎi)=C1ℎi+C2ℎi2+O(ℎi3)其中，C1=u′(x i)，C2=12u″(x i)。
二阶常微分方程的数值求解省公开课获奖课件说课比赛一等奖课件
一. 教学要求
掌握利用降阶把二阶常微分方程转化为一阶微分 方程组，再利用Euler措施数值求解,并能利用MATLAB 软件进行数值计算和符号运算。
二. 教学过程
考虑如下旳二阶微分方程初值问题
d2y dx2 f ( x, y, y ') , y(a) y0 , y '(a) y1, x [a, b]
数值解与真解如下图
例2：利用4阶R-K措施求解例1，并与Euler措施 进行比较。
解 当 h=0.1，即 n=20 时，R-K措施旳Matlab 源程序见
RK_sys1.m，数值成果见下图
rightf_sys1.m
function w=rightf_sys1(x,y,z) w=y;
RK_sys1.m
zk 1
zk
h 6
( L1
2L2
2L3
L4
),
K1 zk , L1 f ( xk , yk , zk ),
K
2
zk
h 2
L1 ,
hh
h
L2 f ( xk 2 , yk 2 K1, zk 2 L1),
K3
zk
h 2
xk 2 , yk 2 K2 , zk 2 L2 ),
>> syms x y >> S=dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0','x')
则能够得到如下旳成果
S=
4/3*cos(x)-1/3*cos(2*x)
注意：只有极少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值措施求数值解。
二阶常微分方程边值问题数值方法
其中 p( x),q( x)为,r已( x知) 函数，则由常微分方程的理论知，通过
变量替换总可以消去方程中的 项，不妨y设 变换后的方程为
y( x) q( x) y( x) r( x)
y(a) ,
y(b)
则近似差分方程成离散差分方程为
yi 1
2 yi h2
yi 1
qi
yi
ri
其中 qi q( xi ), ri r( xi ), i 1,2, , n. y0 ,
第一边界问题：
y0 , yn1
(8.9)
第二边界问题：
y1 y0 h , yn1 yn h
(8.10)
第三边界问题：
y1 (1 0h) y0 1h,
(1 0h) yn1 yn 1h
(8.11)
若 f ( x, y,是y) 的y线, y性 函数时，f 可写成
f (x, y, y) p(x) y( x) q( x) y(x) r( x)
以
y
为待定参数。
0
对第三类边界问题，仍可转化为考虑初值问题(8.5)，取
y0 ,
y0 1 0 y0 ，以 y为0 待定参数。
8.2 有限差分法
将区间[a,b]进行等分：
h
ba, n1
xi
a ih, i 0,1,
,n 1,
设在
x xi , i 0,1, , n 1处的数值解为 。 yi 用中心差分近似微分，即
而且还有误差估
计：
Ri
y( xi )
yi
M 24
h2
(
xi
a)(b xi )
其中 M max y(4。) ( x)
x[a ,b]
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
利用Euler方法求解上述方程组可得如下数 值格式
zk y(a ) y0 , y '(a ) z0 yk yk 1 hzk , zk 1 hf ( xk , yk , zk ), k 1, 2 xk xk 1 h.
其中yk 是y( xk )的近似，zk 是y '( xk )的近似
下的特解，并画出解函数的图形。
>> y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') >> ezplot(y);
例6
d2 y 2 cos( 2 x ) y 求二阶常微分方程 dx 的通解 y(0) 1, y '(0) 0
若令
z y'
, 则上述初值问题可以转化为如下一
阶微分方程组初值问题
x [a , b] y '( x ) z ( x ), z '( x ) f ( x , y( x ), z( x )), x [a , b] z (a ) y z , y(a ) y 1 0 0
当 h=0.1，即 n=20 时，Matlab 源程序见 RK_sys2.m, 数值结
果如下图
rightf_sys1.m
function w=rightf_sys2(x,y,z) w=-y+2*exp(-x)*(x-1);
RK_sys2.m clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=1; Euler_z(1)=1; RK_y(1)=1; RK_z(1)=1; for i=1:length(x)-1 %**** Euler Method ****% Euler_y(i+1)=Euler_y(i)+h*Euler_z(i); Euler_z(i+1)=Euler_z(i)+h*rightf_sys2(x(i),Euler_y(i),Euler_z(i)); %***** R-K4 Method*****% K1=RK_z(i); L1=rightf_sys2(x(i),RK_y(i),RK_z(i)); % K1 and L1
作业 利用Euler方法和R-K方法求解一个 二阶常微分初值问题，并比较数 值结果，计算数值解和解析解的 误差。 利用dsolve函数求解一些微分方程 的通解
数值解与真解如下图
例2：利用4阶R-K方法求解例1，并与Euler方法 进行比较。
解 当 h=0.1，即 n=20 时，R-K方法的Matlab 源程序见
RK_sys1.m，数值结果见下图
rightf_sys1.m
function w=rightf_sys1(x,y,z) w=y;
RK_sys1.m clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=1; Euler_z(1)=-1; %初值 RK_y(1)=1; RK_z(1)=-1; %初值 for i=1:length(x)-1 %**** Euler Method ****% Euler_y(i+1)=Euler_y(i)+h*Euler_z(i); Euler_z(i+1)=Euler_z(i)+h*Euler_y(i);
二阶常微分方程的数值求解
一. 教学要求
掌握利用降阶把二阶常微分方程转化为一阶微分 方程组，再利用Euler方法数值求解,并能利用MATLAB 软件进行数值计算和符号运算。
二. 教学过程
考虑如下的二阶微分方程初值问题
d2 y f ( x , y , y ') , y(a ) y0 , y '(a ) y1 , x [a , b] 2 dx
利用四阶R-K方法求解上述方程组可得如下 数值格式
h yk 1 yk ( K1 2 K 2 2 K 3 K 4 ), 6 h zk 1 zk ( L1 2 L2 2 L3 L4 ), 6 K1 zk , L1 f ( xk , yk , zk ), K 2 zk h L1 , L2 f ( xk h , yk h K1 , zk h L1 ), 2 2 2 2 h h h h K 3 zk L2 , L3 f ( xk , yk K 2 , zk L2 ), 2 2 2 2 K 4 zk hL3 , L4 f ( xk h, yk hK 3 , z k hL3 ).
利用dsolve 函数求微分方程解析解
dsolve 的调用格式
y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v') 其中 y 为输出， eq1、eq2、...为微分方程，cond1、 cond2、...为初值条件，v 为自变量，如果不指定v作为自变 量，则默认t为自变量。
如果省略初值条件，则表示求通解； 如果省略自变量，则默认自变量为
dsolve('Dy=2*x','x'); dsolve('Dy=2*x');
t
％ dy/dx = 2x ％ dy/dt = 2x
若找不到解析解，则返回其积分形式。
例 5：求微分方程 xy ' y e x 0 在初值条件 y(1) 2e
%***** R-K4 Method*****% K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); K2=RK_z(i)+0.5*h*L1;
% K1 and 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 L2 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); % K3 and L3 K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4 RK_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); end plot(x,Euler_y,'r+',x,exp(-x),'k-',x,RK_y,'b*'); xlabel('Variable x'); ylabel('Variable y');
dy x2 2 xy xe 例 4：求微分方程 的通解，并验证。 dx
>> syms x; diff(y)+2*x*y - x*exp(-x^2) >> y=dsolve('Dy+2*x*y=x*exp(-x^2)','x')
几点说明
微分方程中用 D 表示对 自变量 的导数，如： Dy y'； D2y y''； D3y y'''
K2=RK_z(i)+0.5*h*L1; L2=rightf_sys2(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1); % K2 and L2 K3=RK_z(i)+0.5*h*L2; L3=rightf_sys2(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2); % K3 and L3 K4=RK_z(i)+h*L3; L4=rightf_sys2(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4 RK_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); end plot(x,Euler_y,'r+',x,cos(x)+x.*exp(-x),'k-',x,RK_y,'b*'); xlabel('Variable x'); ylabel('Variable y');
clc;clear; h=0.1; a=0;b=2; x=a:h:b; y(1)=1; z(1)=-1; for i=1:length(x)-1 y(i+1)=y(i)+h*z(i); z(i+1)=z(i)+h*y(i); end plot(x,y,'r+',x,exp(-x),'k-'); xlabel('Variable x'); ylabel('Variable y');
x [0, 2] y '( x ) z ( x ), x [0, 2] z '( x ) y( x ), z (0) 1 z , y(0) 1. 0