精华内容
下载资源
问答
  • 常微分方程求解定义
    千次阅读
    2020-12-06 12:12:19

    在惯性导航以及VIO等实际问题中利用IMU求解位姿需要对IMU测量值进行积分得到需要的位置和姿态,其中主要就是求解微分方程。但之前求解微分方程的解析方法主要是应用于一些简单和特殊的微分方程求解中,对于一般形式的微分方程,一般很难用解析方法求出精确解,只能用数值方法求解。该系列主要介绍一些常用的常微分方程的数值解法,主要包括:

    这个系列后面文章会用到前面文章的理论和技术,所以建议按照顺序查看。

    简介

    而在具体求解微分方程中,一般来说给定的条件是
    { x ′ ( t ) = f ( t , x ( t ) ) , a ≤ t ≤ b x ( t 0 ) = x 0 \begin{cases} \mathbf{x}^{\prime}(t)=f(t, \mathbf{x}(t)), \quad a \leq t \leq b \\ \mathbf{x}({t_0}) = x_0 \end{cases} {x(t)=f(t,x(t)),atbx(t0)=x0
    即给定微分方程以及原方程在初始点的值,求原方程在某个 t t t下的 x ( t ) x(t) x(t)原方程的值,这类问题就是(一阶)常微分方程初值求解问题。
    (这里不对常微分方程或者偏微分方程等概念,以及求解微分方程的其他条件如边界条件情况做详细介绍,需要了解的话可以自己google,不影响本系列介绍。)

    解法介绍

    一般来说,(一阶)常微分方程数值解法基本思想是:
    在区间 [ a , b ] [a, b] [a,b]中插入一系列间隔相同为 Δ t \Delta t Δt的离散点
    a ≤ t 0 < t 1 < ⋯ < t i < ⋯ < t n − 1 < t n ≤ b a \leq t_0 < t_1< \cdots < t_i < \cdots < t_{n-1} < t_n \leq b at0<t1<<ti<<tn1<tnb
    常微分方程的数值解法的目的就是在给定合适初始值前提下,建立求解 x ( t ) x(t) x(t)的近似值 x t x_t xt的递推方式,这样求得 x ( t ) x(t) x(t)在各个离散点的近似值。
    给定一个具体的时间间隔 Δ t \Delta_t Δt,可以把该方程转换成差分方程
    x ( t + Δ t ) = x ( t ) + ∫ τ ∈ [ t , t + Δ t ] f ( τ , x ( τ ) ) d τ , a ≤ t ≤ b , (1) \mathbf{x}(t+\Delta t)=\mathbf{x}(t)+\int_{\tau \in [t, t+\Delta t]} f(\tau, \mathbf{x}(\tau)) d \tau, \quad a \leq t \leq b, \tag{1} x(t+Δt)=x(t)+τ[t,t+Δt]f(τ,x(τ))dτ,atb,(1)
    这样利用之前说的离散点可以把 t t t离散化,那么 t i = a + i Δ t t_{i}=a+i \Delta t ti=a+iΔt以及 x i ≈ x ( t i ) \mathbf{x}_{i} \approx \mathbf{x}(t_{i}) xix(ti),所以上面公式 ( 2 ) (2) (2)也可以写成
    x ( t i + 1 ) ≈ x i + 1 = x i + ∫ τ ∈ [ a + i Δ t , a + ( i + 1 ) Δ t ] f ( τ , x ( τ ) ) d τ , (2) \mathbf{x}(t_{i+1}) \approx x_{i+1}=\mathbf{x}_{i}+\int_{\tau \in [a+i \Delta t, a+(i+1) \Delta t]} f(\tau, \mathbf{x}(\tau)) d \tau, \tag{2} x(ti+1)xi+1=xi+τ[a+iΔt,a+(i+1)Δt]f(τ,x(τ))dτ,(2)

    截断误差

    上节介绍常微分方程的数值解法就是利用初值和离散点获得近似值 x t x_t xt,但是和真实值 x ( t ) x(t) x(t)之前还是存在误差,即
    e i = x i − x ( t i ) , (3) \boldsymbol{e}_i = x_i - \mathbf{x}(t_i), \tag{3} ei=xix(ti),(3)
    e i \boldsymbol{e}_i ei则表示该数值解法在离散点 t i t_i ti处的局部截断误差。局部截断误差在一定程度上反应了该数值解法的精度。
    一般来说常用泰勒展开式来计算讨论局部截断误差,后面具体方法会给出具体的对应的截断误差的讨论。这里简单给出定义:
    如果数值解法的局部截断误差为 e i = O ( Δ t P + 1 ) \boldsymbol{e}_i = \boldsymbol{O}(\Delta t^{P+1}) ei=O(ΔtP+1),则称该解法具有P阶精度或P阶解法

    这里 O ( Δ t P + 1 ) \boldsymbol{O}(\Delta t^{P+1}) O(ΔtP+1)为泰勒展开式的余项

    更多相关内容
  • % 模拟一个球从 100 米的高度落下后从楼梯上弹起 米到楼梯上,从 5 米开始。 让楼梯间成为% 由函数 y(x) = 5-floor(x) 定义。 球滚下边缘% 高度为 10 ... % 求解微分方程: [t,y] = SolveOdeWithEvent(odefun, initcon
  • matlab优化常微分方程代码关于这个仓库 这个简单的MATLAB代码是使用四阶Runge-Kutta方法对一阶常微分方程dy / dx = func(x,y)进行数值求解的方法。 由于其简单性,您可以轻松地对其进行修改或将其与其他代码组合...
  • 1.了解积分变换的定义、目的以及意义,熟悉积分变换在求解常微分方程过程中的思想和步骤。 2.总结傅氏变换的基本概念和性质,根据傅氏变换、傅氏逆变换公式和性质举例说明其在求解常微分方程的具体应用,总结归纳出...
  • 耦合常微分方程的数值解。 I&EC 基础。 7:3,1968 年。第 514-517 页 测试问题是一个使用 Stefan-Maxwell 方程和两个 Dirichlet 边界条件的三元扩散问题,来自 Ross Taylor & R. Krishna 撰写的“多组分传质”中的第...
  • 常微分方程及其matlab求解 毕业论文设计 常微分方程及其matlab求解 目 录 摘要1 关键字1 引言1 第一章 一阶微分方程的初等解法1 1.1 变量分离微分方程与变量代换1 1.1.1变量分离微分方程2 1.1.2 可化为变量分离微分...

    41528d3028836879cd698677c3999917.gif常微分方程及其matlab求解 毕业论文设计

    常微分方程及其matlab求解 目 录 摘要1 关键字1 引言1 第一章 一阶微分方程的初等解法1 1.1 变量分离微分方程与变量代换1 1.1.1变量分离微分方程2 1.1.2 可化为变量分离微分方程的类型2 1.2线性分式方程3 1.3线性方程与伯努利方程4 1.4全微分方程与积分因子6 1.4.1 全微分方程6 1.4.2 积分因子7 1.5一阶隐方程8 1.5.1可解出y或x的方程的解法8 1.5.2 不显含y(或x)的方程10 1.6 matlab 在一阶方程中的应用10 第二章 高阶微分方程的解法13 2.1 线性微分方程的一般理论13 1 齐次线性微分方程的解的性质与结构13 2.非奇次线性微分方程与常数变易法14 2.2常系数线性微分方程的解法14 1.复值函数与复值解15 2.3.非齐次线性微分方程的比较系数法18 2.4用matlab解高阶线性方程20 第三章 线性微分方程组24 3.1矩阵指数24 3.2基解矩阵的计算25 1.矩阵A存在n个线性无关的特征向量的情形25 3.3用MATLAB解线性方程组29 致 谢30 参考文献30 Ordinary Differential Equation And Carry On By MATLAB30 Abstract30 Keywords30 常微分方程及其matlab求解 摘要: 本文主要讨论了一阶常微分方程和高阶常微分方程的相关解法问题,并用matlab求解相关方程.文章首先探讨了一阶常微分方程的解法,讨论的主要类型有:变量可分离方程、可化为变量可分离方程的类型、齐次方程、一阶线性微分方程;在解决这些类型的一阶常微分方程时,用到的方法有:变量分离法和一阶线性方程的常数变易法.然后讨论了高阶常微分方程的解法的问题,所讨论的解法有:非齐线性方程的常数变易法、常系数齐线性方程的欧拉待定指数法和非齐线性方程的比较系数法.最后简单讨论了线性微分方程组的解法和性质。 关键字: 一阶常微分方程 高阶常微分方程 解法 matlab Ordinary Differential Equation And Carry On By MATLAB Abstract: This paper mainly discusses some related solutions of the first-order and higher-order ordinary differential equations and solve the relevant equations with matlab. This paper firstly introduces the basic concept of differential equations. on such a basis, the paper probes into the solutions of the first-order differential equations including the main types such as variable separable equation, separable variable equations which can be translated into the equation homogeneous equation,and a linear differential equations. To solve such types of first-order differential equation, the s can be used: variable separation and a linear equation of constant change of law.Then,it discusses the solutions of the higher-order differential equation. The solution are non-homogeneous linear equation of constant variation , Euler be determined index of constant coefficient of linear equations, nonhomogeneous linear equations comparison . Finally,Briefly discussed the solution of linear differential equations and their properties. Keywords: first-order ordinary differential equations; higher-order ordinary differential equations; Solution; MATLAB 引言: 常微分方程作为现代数学的一个重要分支,它的产生几乎与微积分是同时代的,经过历史的演变,它已经是各种应用学科和数学理论研究都不可缺少的工具。随着计算机技术的飞速发展,更使它的应用渗透到力学、天文、物理等各个领域。 遗憾的是,绝大多数微分方程定解问题的解不能以实用的解析形式来表示,这就产生了理论与应用的矛盾:一方面,人们建立了大量实用的数学模型,列出了反映客观现象的微分方程;另一方面,人们又无法得到这些方程的准确解以定量地描述客观过程。从20世纪80年代以来,世界各国所开发的数学类科技软件多达几十种,在我国流行的数学软件主要有四种:MATLAB、Mathematica、Maple和MathCAD。其中,MATLAB有着其它几种数学软件无法比拟的优势和适用面,近几年,MATLAB已成为科学工作者首选的数学软件。本文把微分方程求解和MATLAB有机的结合起来,全面介绍了微分方程的求解在MATLAB中的实现,使得让数学基础不深厚的读者同样能轻易利用MATLAB解决较高深的微分方程问题。 第一章 一阶微分方程的初等解法 1.1 变量分离微分方程与变量代换 变量分离微分方程是一种最简单也是最基本的可用初等方法求解的微分方程类型,对一般的微分方程总是设法寻求适当的变量代换,将其化为变量分离微分方程来求解。 1.1.1变量分离微分方程 形如 (1.1) 的方程,称为变量分离微分方程,其中分别是的连续函数。 如果,方程改写为: (1.2) 如果是方程的解,且,则在解的定义域内满足 (1.3

    展开全文
  • 参考:《高等应用数学问题的MATLAB求解(第四版)》 本文涵盖微分方程的转换、刚性微分方程、隐式微分方程、微分代数方程求解方法。

    目录

    微分方程的转换

    一、单个高阶常微分方程

    二、高阶常微分方程组

    刚性微分方程求解

    隐式微分方程求解

    微分代数方程求解


     

    微分方程的转换

    根据微分方程求解的标准型,要得到微分方程的数值解,应该先将该方程变换成一阶常微分方程组

    一、单个高阶常微分方程

    1、原状态

    高阶常微分方程一般形式:

    gif.latex?y^{(n)}=f(t,y,y{}',...,y^{(n-1)})

    初值为y(0),y'(0),...,y^(n-1)(0)

    2、变换后

    可以选择一组状态变量x1=y,x2=y',...,xn=y^(n-1),可以将原方程变换为一阶方程组形式:

    • x'1=x2
    • x'2=x3
    • ...
    • x'n=f(t,x1,x2,...,xn)

    初值为x1(0)=y(0),x2(0)=y'(0),...,xn(0)=y^(n-1)(0)

    例1.1 求y''+mu*(y^2-1)*y'+y=0,初值为y(0)=0.2,y'(0)=-0.7

    代码如下:

    %% 单个高阶常微分方程处理方法f=@(t,x,mu) [x(2);-mu*(x(1)^2-1)*x(2)-x(1)];% 带附加参数的匿名函数x0=[-0.2;-0.7];tn=20;mu=1;tic; % tic,toc记录时间[t1,y1]=ode45(f,[0,tn],x0,[],mu);toc;mu=2;tic;[t2,y2]=ode45(f,[0,tn],x0,[],mu);toc;figure(1);plot(t1,y1(:,1),t2,y2(:,1),'--');figure(2);plot(t1,y1(:,2),t2,y2(:,2),'--');figure(3);plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),'--'); %绘制相平面图% mu=1000,tn=3000,x0=[2;0]时,未出现内存不够的情况,但用求解刚性微分方程的方式计算量小得多% options=odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);是解微分方程时的选项设置% 'RelTol',1e-4,是相对误差设置,'AbsTol',[1e-4 1e-4 1e-5]是绝对误差设置f=@(t,x,mu) [x(2);-mu*(x(1)^2-1)*x(2)-x(1)];x0=[2;0];tn=3000;mu=1000;tic;%options=odeset('RelTol',1e-10,'AbsTol',1e-10);[t3,y3]=ode45(f,[0,tn],x0,[],mu);toc;figure(4);plot(t3,y3(:,1));figure(5);plot(t3,y3(:,2),'--');figure(6);plot(y3(:,1),y3(:,2)); %绘制相平面图

    运行结果:

    mu=1和mu=2时,实线为[t,y];虚线为[t,y'];用时分别为 0.002099 秒和0.001892 秒。

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

      mu=1000,tn=3000,x0=[2;0]时,实线为[t,y];虚线为[t,y'];用时为28.018461 秒。

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

    注:用ode45()计算mu=1000,tn=3000,x0=[2;0]时的方程,如果电脑内存不够的话会报错,因为步长过小导致,可采取刚性方程的算法,用ode15s可大大减小计算量。

    二、高阶常微分方程组

    1、原状态

    如果可以显式地将两个方程写成:

    • gif.latex?x^{(m)}=f(t,x,x',...,x^{(m-1)},y,y',...y^{(n-1)})
    • gif.latex?y^{(n)}=g(t,x,x',...,x^{(m-1)},y,y',...,y^{(n-1)})

    初值为y(0),y'(0),...,y^(n-1)(0)

    2、变换后

    仍然可以选择一组状态变量x1=x,x2=x',...,xm=x^(m-1),x(m+1)=y,x(m+2)=y',...,x(m+n)=y^(n-1),可以将原方程变换为一阶微分方程组:

    • x'1=x2
    • ...
    • x'm=f(t,x1,x2,...,x(m+n))
    • x'(m+1)=x(m+2)
    • ...
    • x'(m+n)=g(t,x1,x2,...,x(m+n))

    例1.2 卫星的运动轨迹满足下面方程组:

    • x''=2*y'+x-mu1*(x+mu)/r1^3-mu*(x-mu1)/r2^3
    • y''=-2*x'+y-mu1*y/r1^3-mu*y/r2^3

     其中,mu=1/82.45,mu1=1-mu,r1=sqrt((x+mu)^2+y^2),r2=sqrt((x-mu1)^2+y^2)

    代码如下:

    函数程序

    function f=c_7_3_fun% 定义一个接口,通过匿名函数调用% 通过定义接口的方式,在一个函数文件例里可编写多个函数f.apolloeq=@apolloeq;end% 选择一组状态向量:x1=x,x2=x',x3=y,x4=y'function dx=apolloeq(t,x) % 带有附加参数的微分方程描述曲线mu=1/82.45;mu1=1-mu;r1=sqrt((x(1)+mu)^2+x(3)^2);r2=sqrt((x(1)-mu1)^2+x(3)^2);dx=[x(2);2*x(4)+x(1)-mu1*(x(1)+mu)/r1^3-mu*(x(1)-mu1)/r2^3;x(4);-2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];end

    主程序

    %% 高阶常微分方程组变换方法f=c_7_3_funx0=[1.2;0;0;-1.04935751];tic;[t,y]=ode45(@f.apolloeq,[0,20],x0);toc;length(t) % 读取向量长度plot(y(:,1),y(:,3)); % 绘制相平面图,默认精度不够(0.001)%精度改为1e-5tic;[t1,y1]=ode45(@f.apolloeq,[0,20],x0,odeset('RelTol',1e-5));toc;length(t1) % 读取向量长度figure;plot(y1(:,1),y1(:,3)); %精度改为1e-6tic;[t2,y2]=ode45(@f.apolloeq,[0,20],x0,odeset('RelTol',1e-6));toc;length(t2) % 读取向量长度figure;plot(y2(:,1),y2(:,3)); min(diff(t2)) %计算步长figure;plot(t2(1:end-1),diff(t2)); % 绘制实际使用的步长变化曲线%% 采用定步长四阶RK算法f=c_7_3_funx0=[1.2;0;0;-1.04935751];tic;[t3,y3]=rk_4(@f.apolloeq,[0,20,0.01],x0);toc;length(t3) % 读取向量长度figure;plot(y3(:,1),y3(:,3)); % 步长改为0.001f=c_7_3_funx0=[1.2;0;0;-1.04935751];tic;[t4,y4]=rk_4(@f.apolloeq,[0,20,0.001],x0);toc;length(t4) % 读取向量长度figure;plot(y4(:,1),y4(:,3)); %在求解常微分方程时建议采用变步长算法

    运行结果:(算法比较)

    算法:ode45(),默认精度0.001,历时 0.007719 秒,ans =689(向量长度)

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

    算法:ode45(),精度1e-5,历时 0.017587 秒,ans =1309(向量长度)

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

    算法:ode45(),精度1e-6,历时 0.017856 秒,ans =1873(向量长度),ans =1.8927e-04(最小步长)

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

     注:改变算法精度重新求解,发现结果是否不同,可以通过此法验证数值解的正确性。

    算法:定步长四阶RK算法,步长0.01,历时 0.044436 秒,ans =2001(向量长度)结果错误

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

     算法:定步长四阶RK算法,步长0.001,历时 0.716064 秒,ans =20001(向量长度)

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

    注:定步长四阶RK算法中,当步长减小至0.001时,结果才大致正确,大求解时间明显变长,建议采用变步长算法。

    例1.3 将下列方程组转化为一阶微分方程组。

    • x''+2*y'*x=2*y''
    • x''*y'+3*x'*y''+x*y'-y=5 

    代码如下:

    %% 高阶微分方程隐式地含有最高项% 除了手工转化,可以利用MATLAB符号运算工具箱syms x1 x2 x3 x4 dx dy % 声明符号变量[dx,dy]=solve([dx+2*x4*x1==2*dy,dx*x4+3*x2*dy+x1*x4-x3==5],[dx,dy])% 解代数方程% solve(eqn,var)

    运行结果:

    dx =(2*(x3 - x1*x4 - 3*x1*x2*x4 + 5))/(3*x2 + 2*x4)
    dy =(2*x1*x4^2 - x1*x4 + x3 + 5)/(3*x2 + 2*x4)

    注:可以选择手动求解,正常选择状态变量,然后代入法求x''和y'',也可利用solve求代数方程的方法,更高效。 

     

    刚性微分方程求解

    刚性微分方程:一些解变化缓慢,另一些变化快,且相差悬殊。

    注:刚性问题——在用微分方程描述的一个变化过程中,往往又包含着多个相互作用但变化速相差十分悬殊的子过程,这样一类过程就认为具有“刚性”。描述这类过程的微分方程初值问题称为“刚性问题”。例如,宇航飞行器自动控制系统一般包含两个相互作用但效应速度相差十分悬殊的子系统,一个是控制飞行器质心运动的系统,当飞行器速度较大时,质心运动惯性较大,因而相对来说变化缓慢;另一个是控制飞行器运动姿态的系统,由于惯性小,相对来说变化很快,因而整个系统就是一个刚性系统。

    例2.1 (例1.1)求y''+mu*(y^2-1)*y'+y=0,初值为y(0)=0.2,y'(0)=-0.7,用刚性方程解法求解mu=1000,tn=3000,x0=[2;0]的情况。

    常用求解器:ode15s()

    代码如下:

    %% 刚性微分方程求解1% 对于刚性问题的求解,ode15s求解速度比ode45效率高得多ff=odeset;ff.RelTol=1e-10;ff.AbsTol=1e-10; % 设置控制精度x0=[2;0];tn=3000;f=@(t,x,mu) [x(2);-mu*(x(1)^2-1)*x(2)-x(1)];mu=1000;tic;[t4,y4]=ode15s(f,[0,tn],x0,ff,mu);toc;length(t4)plot(t4,y4(:,1));figure;plot(t4,y4(:,2),'--');

    运行结果:(mu=1000,tn=3000,x0=[2;0]时)

    算法:ode15s,实线为[t,y];虚线为[t,y'];用时为0.064681 秒,取点5917个。

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

    对比组(代码见例1.1)

    算法:ode45,实线为[t,y];虚线为[t,y'];用时为28.018461 秒,取点6757169个。

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

     注:1、可见,在对刚性问题的求解中,用专用的刚性求解器(此处为ode15s)能大大提高效率。2、并不是所有刚性问题求解都要刻意选择钢性问题的解法,许多传统的刚性问题采用普通求解方法就可以解出。

    例2.2 考虑下列常微分方程

    • y'1=0.04*(1-y1)-(1-y1)*y1+0.0001*(1-y2)^2
    • y'2=-10^4*y1+3000*(1-y2)^2

    其中,初值为y1(0)=0,y2(0)=1,计算区间为(0,100),选择适当的方法求数值解。

    代码如下:

    %% 刚性微分方程求解2% 用ode45求解f=@(t,y)[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;...    -10^4*y(1)+3000*(1-y(2))^2]; % 用匿名函数描述微分微分方程tic;[t2,y2]=ode45(f,[0,100],[0;1]);toc;length(t2)figure;plot(t2,y2)% 改用ode15s求解f=@(t,y)[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;...    -10^4*y(1)+3000*(1-y(2))^2]; % 用匿名函数描述微分微分方程tic;[t2,y2]=ode15s(f,[0,100],[0;1]);toc;length(t2)figure;plot(t2,y2)% 改用ode15s求解,提高精度ff=odeset;ff.RelTol=1e-10;ff.AbsTol=1e-10; % 设置控制精度f=@(t,y)[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;...    -10^4*y(1)+3000*(1-y(2))^2]; % 用匿名函数描述微分微分方程tic;[t2,y2]=ode15s(f,[0,100],[0;1],ff);toc;length(t2)figure;plot(t2,y2)

    运行结果:

    分别采用ode45、ode15s、提高精度后的ode15s求解。

    ode45(默认精度0.001):历时 0.486060 秒,ans = 356941(点数)

    ode15s(默认精度0.001):历时 0.003206 秒,ans =56(点数)

    ode15s(1e-10):历时 0.012131 秒,ans =614(点数)

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

     

    隐式微分方程求解

    隐式微分方程:不能转换成一阶显式微分方程的微分方程。

    常用求解器:ode15i(),ode15s()(代数方程)

    解法:

    1、矩阵法

    例3.1 给定隐式微分方程组,已知x1(0)=x2(0)=0,求数值解。

    • x'1*sin(x1)+x'2*cos(x2)+x1=1
    • -x'1*cos(x2)+x'2*sin(x1)+x2=0

    解:将其改写成矩阵形式A(x)*x'=B(x)

    • A(x)=[sin(x(1)) cos(x(2));-cos(x(2)) sin(x(1))]
    • B(x)=[1-x(1);-x(2)]

    假设A(x)为非奇异矩阵,则x'=A(X)^(-1)*B(x)(该方法只适合A(x)为非奇异矩阵的情况)

    代码如下:

    %% 隐式微分方程求解% 若得不到奇异矩阵的报错,则正确f=@(t,x)inv([sin(x(1)) cos(x(2));-cos(x(2)) sin(x(1))])*[1-x(1);-x(2)];% inv矩阵求逆opt=odeset;opt.RelTol=1e-6;[t,x]=ode45(f,[0,10],[0;0],opt);plot(t,x)% 此解法有隐患,参照代数方程解法可以更好的求解

    运行结果:

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

    注:此解法有隐患(需要非奇异矩阵的条件),参照下一章 “代数方程解法” 可以更好的求解(例4.2)

    2、巧用代数方程

    例3.2 考虑下列复杂隐式方程

    • x''sin(y')+y''^2=-2*x*y*e^(-x')+x*x''*y'
    • x*x''*y''+cos(y'')=3*y*x'*e^(-x)

    选定状态变量x1=x,x2=x',x3=y,x4=y',设p1=x'',p2=y'',原方程组变为:

    • p(1)*sin(x(4))+p(2)^2+2*x(1)*x(3)*exp(-x(2))-x(1)*p(1)*x(4)=0
    • x(1)*p(1)*p(2)+cos(p(2))-3*x(3)*x(2)*exp(-x(1)) =0

     代码如下:

    函数程序

    function f=c_7_3_fun% 定义一个接口,通过匿名函数调用% 通过定义接口的方式,在一个函数文件例里可编写多个函数%f.apolloeq=@apolloeq;f.c7impode=@c7impode;endfunction dy=c7impode(t,x) % x为参数dx=@(p,x) [p(1)*sin(x(4))+p(2)^2+2*x(1)*x(3)*exp(-x(2))-x(1)*...    p(1)*x(4);x(1)*p(1)*p(2)+cos(p(2))-3*x(3)*x(2)*exp(-x(1))];    % 用x1,x2,x3,x4表示p1,p2ff=optimset;ff.Display='off';% Level of display (see Iterative Display):    % 'off' or 'none' displays no output.    % 'iter' displays output at each iteration,    %        and gives the default exit message.    % 'iter-detailed' displays output at each iteration,     %        and gives the technical exit message.    % 'final' (default) displays just the final output,     %        and gives the default exit message.    % 'final-detailed' displays just the final output,    %        and gives the technical exit message.dx1=fsolve(dx,x([1,3]),ff,x);    % x = fsolve(fun,x0,options)    % fsolve求解非线性方程,初值选择p1(0)=x1,p2(0)=x3,加快代数方程收敛速度和精度    % x作为已知附加参数给出,得出结果为p1,p2dy=[x(2);dx1(1);x(4);dx1(2)];    % dy=[x2=x',dx1(1)=p1=x'',x4=y',dx1(2)=p2=y'']    % 求ode对应的解为[x,x',y,y']end

    主程序

    %% 复杂隐式微分方程f=c_7_3_fun[t,x]=ode15s(f.c7impode,[0,2],[1,0,0,1]);% x1随t变化的图像figure(1);plot(t,x(:,1));% x2随t变化的图像figure(2);plot(t,x(:,2));% x3随t变化的图像figure(3);plot(t,x(:,3));% x4随t变化的图像figure(4);plot(t,x(:,4));% x随t变化的总图像figure(5);plot(t,x);% 下面代码找不到显式解(方程和系统求解器无解)% syms x1 x2 x3 x4 dx dy % 声明符号变量% [dx,dy]=solve([dx*sin(x4)+dy^2+2*x1*x3*exp(-x2)-x1*...%     dx*x4==0,x1*dx*dy+cos(dy)-3*x3*x2*exp(-x1)==0],[dx,dy])

     运行结果:

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

     watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

    3、直接用ode15i()函数求解

     求解之前需要给出(x0,x'0),不能随意赋值,只能有n个是独立的,其余的需要用隐式方程求解,否则将出现矛盾的初始条件,用decic()函数能的处相容的初值。

    例 3.3 (例3.2)考虑下列复杂隐式方程

    • x''sin(y')+y''^2=-2*x*y*e^(-x')+x*x''*y'
    • x*x''*y''+cos(y'')=3*y*x'*e^(-x)

    选定状态变量x1=x,x2=x',x3=y,x4=y',原方程组变为隐式方程标准型:

    • x'1-x2=0
    • x'2*sin(x4)+x'4^2+2*e^(-x2)*x1*x3-x1*x'2*x4=0
    • x'3-x4=0
    • x1*x'2*x'4+cosx'4-3*e^(-x1)*x3*x2=0

    代码如下:

    %% 利用ode15i()求解复杂隐式微分方程% 隐式微分方程的描述f=@(t,x,xd)[xd(1)-x(2);            xd(2)*sin(x(4))+xd(4)^2+2*exp(-x(2))*x(1)*x(3)-x(1)*xd(2)*x(4);            xd(3)-x(4);            x(1)*xd(2)*xd(4)+cos(xd(4))-3*exp(-x(1))*x(3)*x(2)];x0=[1,0,0,1]; % 初值xd0=[1,2,3,-5]; % 无需给出确定值,但要判断正负x0f=[1 1 1 1]; % 保留x0xd0f=[]; % 重新计算xd[x0,xd0]=decic(f,0,x0,x0f,xd0,xd0f)% decic为ode15i计算一致的初始条件% t0=0% 求出相容的初值,由x0确定x0'[t,x]=ode15i(f,[0,2],x0,xd0);plot(t,x) % 绘制时间响应曲线

     运行结果:(同例3.2)

    x0 =       1     0     0     1
    xd0 =  -0.0000  1.6833    1.0000   -0.5166

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

     

    微分代数方程求解

    微分代数方程:指在微分方程中,某些变量间满足满足某些代数方程的约束。

    常用求解器:ode15s()、ode45()

    微分方程更一般的形式可以写成

    M(t,x)*x'=f(t,x)

    对于真正的微分代数方程,M(t,x)为奇异矩阵

    例4.1 考虑下面微分代数方程

    • x'1=-0.2*x1+x2*x3+0.3*x1*x2
    • x'2=2*x1*x2-5*x2*x3-2*x2*x2
    • 0=x1+x2+x3-1

     初始条件:x1(0)=0.8,x2(0)=x3(0)=0.1

    解:用矩阵形式表示该微分代数方程

    [1 0 0;0 1 0;0 0 0]*[x'1 x'2 x'3]=[-0.2*x1+x2*x3+0.3*x1*x2;2*x1*x2-5*x2*x3-2*x2*x2;x1+x2+x3-1]

    代码如下:

    %% 微分代数方程求解f=@(t,x)[-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);    2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);    x(1)+x(2)+x(3)-1];M=[1,0,0;0,1,0;0,0,0]; % 输入质量矩阵options=odeset;options.Mass=M; % 设置质量矩阵x0=[0.8;0.1;0.1];[t,x]=ode15s(f,[0,20],x0,options);plot(t,x)% 也可以转换成常微分方程求解fDae=@(t,x)[-0.2*x(1)+x(2)*(1-x(1)-x(2))+0.3*x(1)*x(2);    2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2))-2*x(2)*x(2)]; % 方程右侧的形式模型x0=[0.8,0.1]; % 初值[t,x]=ode45(fDae,[0,20],x0);x3=1-x(:,1)-x(:,2);figure;plot(t,x,t,x3);% 也可以用隐式微分方程求解的方法求解f=@(t,x,xd)[xd(1)+0.2*x(1)-x(2)*x(3)-0.3*x(1)*x(2);    xd(2)-2*x(1)*x(2)+5*x(2)*x(3)+2*x(2)*x(2);    x(1)+x(2)+x(3)-1]; % 隐式微分方程x0=[0.8,0.1,2];x0f=[1,1,0];xd0=[1,1,1];xd0f=[];[x0,xd0]=decic(f,0,x0,x0f,xd0,xd0f); % 相容初始条件% 隐式微分方程求解与绘图res=ode15i(f,[0,20],x0,xd0);% res为struct结构体数组plot(res.x,res.y) 

    运行结果:

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

    例4.2 (例3.1 给定隐式微分方程组,已知x1(0)=x2(0)=0,求数值解。

    • x'1*sin(x1)+x'2*cos(x2)+x1=1
    • -x'1*cos(x2)+x'2*sin(x1)+x2=0

     代码如下:

    f=@(t,x)[1-x(1);-x(2)];M=@(t,x)[sin(x(1)),cos(x(2));-cos(x(2)),sin(x(1))]; % 质量矩阵options=odeset;options.Mass=M;options.RelTol=1e-6; % 设置求解精度[t,x]=ode45(f,[0,10],[0;0],options);plot(t,x)

    运行结果:

    watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBA5rex5rW35rex5aSc5rex,size_20,color_FFFFFF,t_70,g_se,x_16

     

     

    展开全文
  • 高阶线性常微分方程求解

    3.1 线性常微分方程通解的结构

    n阶线性常微分方程的一般形式为:
    y ( k ) + p 1 ( x ) y ( n − 1 ) + . . . + p n − 1 ( x ) y ′ + p n ( x ) y = f ( x ) , ( 3 − 1 ) y^{(k) }+ p_1(x)y^{(n-1)}+...+ p_{n-1}(x)y^{'}+ p_n(x)y = f(x), \qquad (3-1) y(k)+p1(x)y(n1)+...+pn1(x)y+pn(x)y=f(x),(31) y ( k ) + p 1 ( x ) y ( n − 1 ) + . . . + p n − 1 ( x ) y ′ + p n ( x ) y = 0 ( 齐 次 ) , ( 3 − 2 ) y^{(k) }+ p_1(x)y^{(n-1)}+...+ p_{n-1}(x)y^{'}+ p_n(x)y = 0(齐次), \qquad (3-2) y(k)+p1(x)y(n1)+...+pn1(x)y+pn(x)y=0,(32)

    3.1.1. 齐次高阶通解(讨论了解的性质,并未求解)

    定理1: 叠加原理
    如果 y 1 ( x ) , y 2 ( x ) , . . . , y k ( x ) y_1(x),y_2(x),...,y_k(x) y1(x),y2(x),...,yk(x)为方程(3-2)的k个解,则对任意常数: C 1 , C 2 , . . . , C k C_1,C_2,...,C_k C1,C2,...,Ck,函数
    y = C 1 y 1 ( x ) + C 2 y 2 ( x ) + . . . + C k y k ( x ) y= C_1y_1(x)+ C_2y_2(x)+...+C_ky_k(x) y=C1y1(x)+C2y2(x)+...+Ckyk(x)也是方程的解。

    定义1: y 1 ( x ) , y 2 ( x ) , . . . , y k ( x ) y_1(x),y_2(x),...,y_k(x) y1(x),y2(x),...,yk(x)定义在区间I上,存在不全为零的常数 C 1 , C 2 , . . . , C k C_1,C_2,...,C_k C1,C2,...,Ck使得 ∀ x ∈ I \forall x \in I xI
    C 1 y 1 ( x ) + C 2 y 2 ( x ) + . . . + C k y k ( x ) = 0 C_1y_1(x)+ C_2y_2(x)+...+C_ky_k(x)= 0 C1y1(x)+C2y2(x)+...+Ckyk(x=0则称这n个函数在区间I上是线性相关的,否则线性无关。

    定义2: 由定义在区间I上的n个n-1次函数 y 1 ( x ) , y 2 ( x ) , … , y n ( x ) y_1(x),y_2(x),\ldots,y_n(x) y1(x),y2(x),,yn(x)构成的行列式
    W ( x ) = ∣ y 1 ( x ) y 2 ( x ) ⋯ y k ( x ) y 1 ′ ( x ) y 2 ′ ( x ) ⋯ y k ′ ( x ) ⋮ ⋮ ⋯ ⋮ y 1 n − 1 ( x ) y 2 n − 1 ( x ) ⋯ y k n − 1 ( x ) ∣ W(x)=\left | \begin{array}{cccc} y_1(x) & y_2(x) & \cdots & y_k(x)\\ y_1'(x) & y_2'(x) & \cdots & y_k'(x)\\ \vdots & \vdots & \cdots & \vdots\\ y_1^{n-1}(x) & y_2^{n-1}(x) & \cdots & y_k^{n-1}(x)\\ \end{array} \right| W(x)=y1(x)y1(x)y1n1(x)y2(x)y2(x)y2n1(x)yk(x)yk(x)ykn1(x)称为这n个函数的朗斯基行列式(Wronsky).

    定理2: n阶齐次线性方程(3-2)的n个解 y 1 ( x ) , y 2 ( x ) , … , y k ( x ) y_1(x),y_2(x),\ldots,y_k(x) y1(x),y2(x),,yk(x)定义在区间I上线性无关的充要条件是在I上存在 x 0 x_0 x0 使得他们的朗斯基行列式 W ( x 0 ) ≠ 0 W(x_0) \neq 0 W(x0)=0 .

    定理3: 若n个解 y 1 ( x ) , y 2 ( x ) , … , y k ( x ) y_1(x),y_2(x),\ldots,y_k(x) y1(x),y2(x),,yk(x)是方程(3-2)在区间I上的线性无关解,那么(3-2)的通解可表示为:
    y = C 1 y 1 ( x ) + C 2 y 2 ( x ) + . . . + C n y k ( n ) , ( 3 − 3 ) y= C_1y_1(x)+ C_2y_2(x)+...+C_ny_k(n),\quad (3-3) y=C1y1(x)+C2y2(x)+...+Cnyk(n),(33)它包含了所有解。

    推论: 方程(3-2)线性无关解组中包含的解的个数最大为n 个。

    这n个线性无关解称为方程的基本解组。

    例: y ′ ′ + p 1 ( x ) y ′ + p 2 y = 0 , 已 知 其 一 个 解 为 y 1 , 求 另 一 个 线 性 无 关 解 . y'' + p_1(x)y' + p_2y = 0\quad ,已知其一个解为y_1,求另一个线性无关解. y+p1(x)y+p2y=0y1线.
    设 : y 2 = C ( x ) y 1 , 设:y_2 = C(x)y_1, \\ y2=C(x)y1, y 2 ′ = C ( x ) ′ y 1 + C ( x ) y 1 ′ , y_2' = C(x)'y_1 + C(x)y_1',\\ y2=C(x)y1+C(x)y1, y 2 ′ ′ = C ( x ) ′ ′ y 1 + 2 C ( x ) ′ y 1 ′ + C ( x ) y 1 ′ ′ . y_2'' = C(x)''y_1 + 2C(x)'y_1' + C(x)y_1''. y2=C(x)y1+2C(x)y1+C(x)y1. 带 入 方 程 得 : y 1 C ( x ) ′ ′ + ( 2 y 1 ′ + p 1 y 1 ) C ( x ) ′ + ( y 1 ′ ′ + p 1 y 1 ′ + p 2 ) C ( x ) = 0 , 带入方程得: y_1C(x)'' + (2y_1' + p_1y_1)C(x)' + (y_1'' + p_1y_1' +p_2)C(x) = 0, :y1C(x)+(2y1+p1y1)C(x)+(y1+p1y1+p2)C(x)=0 得 : y 1 C ( x ) ′ ′ + ( 2 y 1 ′ + p 1 y 1 ) C ( x ) ′ = 0 , 变 量 分 离 法 可 得 得:y_1C(x)'' + (2y_1' + p_1y_1)C(x)' = 0,变量分离法可得 y1C(x)+(2y1+p1y1)C(x)=0 y 2 = C ( x ) y 1 = y 1 ∫ e − ∫ p 1 d x y 1 2 d x , ( 3 − 4 ) y_2 = C(x)y_1 = y_1\int \frac{e^{-\int p_1dx}}{y_1^2}dx,\quad(3-4) y2=C(x)y1=y1y12ep1dxdx,(34) 显 然 y 1 , y 2 线 性 无 关 显然y_1,y_2线性无关 \\ y1,y2线

    3.1.2.非齐次高阶通解(介绍了通用的已知齐次基本解组求解特解的方法)

    定理4: 设y*是方程(3-1)的一个特解, y ‾ \overline{y} y是(3_2)的通解,则(3-1)的通解是 y ∗ + y ‾ y^* + \overline{y} y+y.

    \qquad 上一篇文章2.2中我们使用常数变易法求解了一阶非齐次常微分方程的解(在求出齐次解的情况下,推导非齐次);对于高阶常微分方程我们采取相同的策略,由齐次基本解组推导出非齐次的特解。

    (1)求 C i ( x ) C_i(x) Ci(x)

    设方程(3-2)的通解是(3-3),我们假设特解的形式是将常数 C i C_i Ci换为 C i ( x ) C_i(x) Ci(x),即:
    y = C 1 ( x ) y 1 ( x ) + C 2 ( x ) y 2 ( x ) + . . . + C n ( x ) y k ( n ) , ( 3 − 5 ) y= C_1(x)y_1(x)+ C_2(x)y_2(x)+...+C_n(x)y_k(n),\quad (3-5) y=C1(x)y1(x)+C2(x)y2(x)+...+Cn(x)yk(n),(35)
    \qquad 求解出 C i ( x ) C_i(x) Ci(x)就可以确定特解,由于存在n个待定函数,因此求解还需要n-1个限制条件,这些条件按理说可以随意设置(这个地方需要理解一下随意设置的含义是什么:我的理解是就像是对于方程 f ( x ) = x , 我 们 在 实 数 R 范 围 内 可 以 令 x = s h ( t ) , x = t a n ( t ) 等 等 f(x) = x,我们在实数R范围内可以令x = sh(t),x = tan(t)等等 f(x)=xRx=sh(t)x=tan(t)).

    \qquad 当然随意设置之后得到的方程越简单越好。这里采用的方式是对通解多次求导的方式来获得n-1个条件如下,在前n-1次求导的过程中将含有 C i ′ ( x ) C_i'(x) Ci(x)的项求和设为0
    y ′ = C 1 ( x ) y 1 ′ ( x ) + C 2 ( x ) y 2 ′ ( x ) + . . . + C n ( x ) y k ′ ( n ) + 0 , ( 1 ) y'= C_1(x)y_1'(x)+ C_2(x)y_2'(x)+...+C_n(x)y_k'(n) + 0,\quad (1) y=C1(x)y1(x)+C2(x)y2(x)+...+Cn(x)yk(n)+0,(1) 0 = C 1 ′ ( x ) y 1 ( x ) + C 2 ′ ( x ) y 2 ( x ) + . . . + C n ′ ( x ) y k ( n ) , ( 1.1 ) 0 = C_1'(x)y_1(x)+ C_2'(x)y_2(x)+...+C_n'(x)y_k(n),\quad (1.1) 0=C1(x)y1(x)+C2(x)y2(x)+...+Cn(x)yk(n),(1.1) y ′ ′ = C 1 ( x ) y 1 ′ ′ ( x ) + C 2 ( x ) y 2 ′ ′ ( x ) + . . . + C n ( x ) y k ′ ′ ( n ) + 0 , ( 2 ) y''= C_1(x)y_1''(x)+ C_2(x)y_2''(x)+...+C_n(x)y_k''(n) + 0,\quad (2) y=C1(x)y1(x)+C2(x)y2(x)+...+Cn(x)yk(n)+0,(2) 0 = C 1 ′ ( x ) y 1 ′ ( x ) + C 2 ′ ( x ) y 2 ′ ( x ) + . . . + C n ′ ( x ) y k ′ ( n ) , ( 2.1 ) 0 = C_1'(x)y_1'(x)+ C_2'(x)y_2'(x)+...+C_n'(x)y_k'(n),\quad (2.1) 0=C1(x)y1(x)+C2(x)y2(x)+...+Cn(x)yk(n),(2.1) ⋮ \vdots y n − 1 = C 1 y 1 n − 1 ( x ) + C 2 y 2 n − 1 ( x ) + . . . + C n y k n − 1 ( n ) + 0 , ( n − 1 ) y^{n-1}= C_1y_1^{n-1}(x)+ C_2y_2^{n-1}(x)+...+C_ny_k^{n-1}(n) + 0,\quad (n-1) yn1=C1y1n1(x)+C2y2n1(x)+...+Cnykn1(n)+0,(n1) 0 = C 1 ′ ( x ) y 1 ′ ( x ) + C 2 ′ ( x ) y 2 ′ ( x ) + . . . + C n ′ ( x ) y k ′ ( n ) , ( n − 1.1 ) \qquad0 = C_1'(x)y_1'(x)+ C_2'(x)y_2'(x)+...+C_n'(x)y_k'(n),\quad (n-1.1) 0=C1(x)y1(x)+C2(x)y2(x)+...+Cn(x)yk(n),(n1.1) y n = C 1 ( x ) y 1 n ( x ) + C 2 ( x ) y 2 n ( x ) + . . . + C n ( x ) y k n ( n ) + y^{n}= C_1(x)y_1^{n}(x)+ C_2(x)y_2^{n}(x)+...+C_n(x)y_k^{n}(n)+\qquad yn=C1(x)y1n(x)+C2(x)y2n(x)+...+Cn(x)ykn(n)+ C 1 ′ ( x ) y 1 n − 1 ( x ) + C 2 ′ ( x ) y 2 n − 1 ( x ) + . . . + C n ′ ( x ) y k n − 1 ( n ) , ( n ) C_1'(x)y_1^{n-1}(x)+ C_2'(x)y_2^{n-1}(x)+...+C_n'(x)y_k^{n-1}(n),\quad (n) C1(x)y1n1(x)+C2(x)y2n1(x)+...+Cn(x)ykn1(n),(n)将上面式子(i)带入方程(3-1),每个式子的对应位置的项相加求和为(3-2)解,故获得最终化简式:
    C 1 ′ ( x ) y 1 n − 1 ( x ) + C 2 ′ ( x ) y 2 n − 1 ( x ) + . . . + C n ′ ( x ) y k n − 1 ( n ) = f ( x ) , ( n . 1 ) . C_1'(x)y_1^{n-1}(x)+ C_2'(x)y_2^{n-1}(x)+...+C_n'(x)y_k^{n-1}(n) = f(x),\quad (n.1). C1(x)y1n1(x)+C2(x)y2n1(x)+...+Cn(x)ykn1(n)=f(x),(n.1).由(i.1)组成n个微分方程组:
    ∣ 0 0 ⋮ f ( x ) ∣ = ∣ y 1 ( x ) y 2 ( x ) ⋯ y k ( x ) y 1 ′ ( x ) y 2 ′ ( x ) ⋯ y k ′ ( x ) ⋮ ⋮ ⋯ ⋮ y 1 n − 1 ( x ) y 2 n − 1 ( x ) ⋯ y k n − 1 ( x ) ∣ ⋅ ∣ C 1 ′ ( x ) C 2 ′ ( x ) ⋮ C n ′ ( x ) ∣ = A ⋅ B ( 3 − 6 ) \left | \begin{array}{cccc} 0\\ 0\\ \vdots \\ f(x)\\ \end{array} \right| =\left | \begin{array}{cccc} y_1(x) & y_2(x) & \cdots & y_k(x)\\ y_1'(x) & y_2'(x) & \cdots & y_k'(x)\\ \vdots & \vdots & \cdots & \vdots\\ y_1^{n-1}(x) & y_2^{n-1}(x) & \cdots & y_k^{n-1}(x)\\ \end{array} \right| ·\left | \begin{array}{cccc} C_1'(x)\\ C_2'(x)\\ \vdots\\ C_n'(x)\\ \end{array} \right| =A·B\quad(3-6) 00f(x)=y1(x)y1(x)y1n1(x)y2(x)y2(x)y2n1(x)yk(x)yk(x)ykn1(x)C1(x)C2(x)Cn(x)=AB(36)显然A是Wronsky行列式,由于基本解组是线性无关的, ∣ A ∣ ≠ 0 |A| \neq 0 A=0,(3-6)有唯一解。

    C i ′ ( x ) = φ i ( x ) C_i'(x) = \varphi_i(x) Ci(x)=φi(x),积分可得 C i ( x ) C_i(x) Ci(x).

    高阶非齐次常微分方程的解为:
    y = ∑ i = 1 n k i y i ( x ) + ∑ i = 1 n y i ( x ) ∫ φ i ( x ) d x . y = \sum_{i = 1}^nk_iy_i(x) + \sum_{i = 1}^ny_i(x)\int \varphi_i(x)dx. y=i=1nkiyi(x)+i=1nyi(x)φi(x)dx.

    (2)例程

    在这里插入图片描述
    在这里插入图片描述

    3.2 常系数齐次线性常微分方程(由二阶常微分方程求解推广到高阶)

    \qquad 在3.1节中我们介绍了高阶常微分方程的解的结构,但是从上面两个例子来看高阶常微分方程的求解是比较麻烦的。我们先解决二阶常微分方程,并将其解法类似推广到高阶常微分方程。

    3.2.1. 二阶常系数齐次微分方程求解

    \qquad 二阶常微分方程的一般形式为:
    y ′ ′ + p 1 y ′ + p 2 y = 0. ( 3 − 7 ) y'' + p_1y' +p_2y = 0.\qquad (3-7) y+p1y+p2y=0.(37)
    我们设 y = e r x y=e^{rx} y=erx,它的各阶导数只相差一个因子r;将其代入方程(3-7),有:
    y ( r 2 + p 1 ( x ) r + 1 ) = 0 , y ≢ 0 y(r^2 + p_1(x)r + 1) = 0,\quad y \not\equiv 0 y(r2+p1(x)r+1)=0,y0 故 得 到 特 征 方 程 : r 2 + p 1 r + 1 ≡ 0. 故得到特征方程:\quad r^2 + p_1r + 1 \equiv 0. r2+p1r+10.
    (1) r 1 ≠ r 2 , y = C 1 e r 1 x + C 2 e r 2 x r_1 \neq r_2 ,y = C_1e^{r_1x} + C_2e^{r_2x} r1=r2y=C1er1x+C2er2x
    (2) r 1 = r 2 , y = ( C 1 + C 2 x ) e r 1 x r_1 = r_2 ,y = (C_1 + C_2x)e^{r_1x} r1=r2y=(C1+C2x)er1x,利用了3.1中例题的结论。
    (3) r 1 , r 2 r_1,r_2 r1,r2是一对共轭虚根 α ± i β \alpha \pm i\beta α±iβ时,y有(1)的形式,即:
    y 1 = e a x ( c o s β x + i s i n β x ) y 2 = e a x ( c o s β x − i s i n β x ) y_1 = e^{ax}(cos \beta x + isin \beta x)\quad y_2 = e^{ax}(cos \beta x - isin \beta x) y1=eax(cosβx+isinβx)y2=eax(cosβxisinβx) y = C 1 e r 1 x + C 2 e r 2 x y = C_1e^{r_1x} + C_2e^{r_2x} y=C1er1x+C2er2x 令 C 3 = ( C 1 + C 2 ) / 2 , C 4 = ( C 1 − C 2 ) / 2 i 令C_3 = (C_1 + C_2)/2,\quad C_4 = (C_1 - C_2)/2i C3=(C1+C2)/2,C4=(C1C2)/2i y = e a x ( C 3 c o s β x + C 4 s i n β x ) , C i 是 复 数 y = e^{ax}(C_3 cos \beta x + C_4sin \beta x),\quad C_i是复数 y=eax(C3cosβx+C4sinβx),Ci

    3.2.2. 推广到高阶

    高阶时我们同样列出特征方程,分析根的组成结构类似二阶常系数齐次微分方程,我们可以得到:
    在这里插入图片描述

    3.3 常系数非齐次线性常微分方程

    3.3.1 待定系数法(仍然以二阶为例,推广高阶)

    \qquad 二阶常微分方程的一般形式为:
    y ′ ′ + p 1 y ′ + p 2 y = f ( x ) . ( 3 − 8 ) y'' + p_1y' +p_2y = f(x).\qquad (3-8) y+p1y+p2y=f(x).(38)

    类型Ⅰ

    f ( x ) = e λ x P m , P m 是 m 次 多 项 式 ; y ∗ = e λ x Q ( x ) 为 方 程 的 特 解 , Q ( x ) 为 多 项 式 f(x) = e^{\lambda x} P_m,\quad P_m是m次多项式;y^* = e^{\lambda x} Q(x)为方程的特解,Q(x)为多项式 f(x)=eλxPm,Pmmy=eλxQ(x)Q(x)

    f ( x ) f(x) f(x) y ∗ y^* y带入方程(3-8)有:

    ( λ 2 + p 1 λ + p 2 ) Q ( x ) + ( 2 λ + p 1 ) Q ′ ( x ) + Q ′ ′ ( x ) = P m ( x ) , ( 3 − 9 ) (\lambda ^2 + p_1\lambda + p_2)Q(x) + (2\lambda + p_1)Q'(x) + Q''(x) = P_m(x),\quad (3-9) (λ2+p1λ+p2)Q(x)+(2λ+p1)Q(x)+Q(x)=Pm(x),(39)

    • λ \lambda λ 不是方程 λ 2 + p 1 λ + p 2 = 0 \lambda ^2 + p_1\lambda + p_2 = 0 λ2+p1λ+p2=0的根, Q ( x ) Q(x) Q(x) P m ( x ) P_m(x) Pm(x)次数相同,将 λ \lambda λ带入方程(3-9),求取对应幂次的系数即可。
    • λ \lambda λ 是方程的单根时,式(3-9)变为
      ( 2 λ + p 1 ) Q ′ ( x ) + Q ′ ′ ( x ) = P m ( x ) , (2\lambda + p_1)Q'(x) + Q''(x) = P_m(x), (2λ+p1)Q(x)+Q(x)=Pm(x), 设 Q ( x ) = x Q m ( x ) , 求 取 对 应 幂 次 的 系 数 即 可 。 设Q(x) = xQ_m(x),求取对应幂次的系数即可。 Q(x)=xQm(x) 这 里 实 际 上 是 Q ( x ) = ( a x + b ) Q m ( x ) , 但 是 一 次 求 导 之 后 等 效 Q ( x ) = x Q m ( x ) , 即 省 略 了 常 数 项 。 这里实际上是Q(x) = (ax + b)Q_m(x),但是一次求导之后等效Q(x) = xQ_m(x),即省略了常数项。 Q(x)=(ax+b)Qm(x)Q(x)=xQm(x)
    • λ \lambda λ 是方程的重根时,式(3-9)变为
      Q ′ ′ ( x ) = P m ( x ) , Q''(x) = P_m(x), Q(x)=Pm(x), 同 理 设 Q ( x ) = x 2 Q m ( x ) , 求 取 对 应 幂 次 的 系 数 即 可 。 同理设Q(x) = x^2Q_m(x),求取对应幂次的系数即可。 Q(x)=x2Qm(x)

    该结论可以推广到高阶,例如下:
    在这里插入图片描述

    类型Ⅱ

    f ( x ) = e a x [ P l ( x ) c o s β x + Q m ( x ) s i n ( β x ) ] , P l 是 l 次 多 项 式 , Q m 是 m 次 多 项 式 f(x) = e^{ax} [P_l(x)cos \beta x + Q_m(x)sin(\beta x)],\quad P_l是l次多项式,Q_m是m次多项式 f(x)=eax[Pl(x)cosβx+Qm(x)sin(βx)],Pll,Qmm

    1.这种形式的方程有如下特解:
    y ∗ = x ∗ k e a x [ R n ( 1 ) ( x ) c o s β x + R n ( 2 ) ( x ) s i n β x ] , R n ( x ) 是 n 次 多 项 式 y^* = x*ke^{ax}[R_n^{(1)}(x)cos\beta x + R_n^{(2)}(x)sin\beta x],R_n(x)是n次多项式 y=xkeax[Rn(1)(x)cosβx+Rn(2)(x)sinβx]Rn(x)n λ = α ± i β 是 方 程 的 根 的 时 候 k = 1 , 否 则 为 0 ; n = m a x { l , m } \lambda = \alpha \pm i\beta是方程的根的时候k= 1,否则为0;n = max\{l,m\} λ=α±iβk=10n=max{l,m}
    2.书中该结论的证明以复杂为理由省略,我自己粗略证明了一下(思路就是将类型Ⅰ扩充到复数域):

    • 先求解下面方程的解:
      当 f ( x ) = e α x P l ( x ) c o s β x = e ( α + i β ) x + e ( α + i β ) x 2 P l ( x ) , λ 1 = α + i β , λ 2 = α − i β . 当f(x) = e^{\alpha x} P_l(x)cos\beta x =\frac{e^{(\alpha + i\beta)x} + e^{(\alpha + i\beta)x}}2P_l(x),\quad \lambda _1= \alpha + i\beta,\lambda_2 = \alpha - i\beta. f(x)=eαxPl(x)cosβx=2e(α+iβ)x+e(α+iβ)xPl(x),λ1=α+iβ,λ2=αiβ. 根 据 类 型 Ⅰ 中 的 分 析 , 我 们 可 以 知 道 λ 1 , 2 只 可 能 是 方 程 λ 2 + p 1 λ + p 2 = 0 的 单 根 或 者 不 是 根 , 因 为 λ 是 虚 数 。 根据类型Ⅰ中的分析,我们可以知道\lambda _1,_2 只可能是方程\lambda ^2 + p_1\lambda + p_2 = 0的单根或者不是根,因为\lambda 是虚数。 λ1,2λ2+p1λ+p2=0λ 设 特 解 为 : y 1 ∗ = x k ( R l ( 1 ) ( x ) e λ 1 x + ( R l ( 2 ) ( x ) e λ 2 x ) 设特解为:y_1^* = x^k(R_l^{(1)}(x) e^{\lambda _1 x} + (R_l^{(2)}(x) e^{\lambda _2 x}) y1=xk(Rl(1)(x)eλ1x+(Rl(2)(x)eλ2x)

    • 同理当 g ( x ) = e a x Q m ( x ) s i n β x g(x) = e^{ax} Q_m(x)sin\beta x g(x)=eaxQm(x)sinβx,有以下形式的特解:
      y 2 ∗ = x k ( R m ( 1 ) ( x ) e λ 1 x + ( R m ( 2 ) ( x ) e λ 2 x ) . y_2^* = x^k(R_m^{(1)}(x) e^{\lambda _1 x} + (R_m^{(2)}(x) e^{\lambda _2 x}). y2=xk(Rm(1)(x)eλ1x+(Rm(2)(x)eλ2x).

    • 故高阶常微分方程等号右边为f(x)+g(x)时,有特解:
      y ∗ = x k ( R l ( 1 ) ( x ) e λ 1 x + ( R l ( 2 ) ( x ) e λ 2 x + R m ( 1 ) ( x ) e λ 1 x + ( R m ( 2 ) ( x ) e λ 2 x ) . y^* = x^k(R_l^{(1)}(x) e^{\lambda _1 x} + (R_l^{(2)}(x) e^{\lambda _2 x} + R_m^{(1)}(x) e^{\lambda _1 x} + (R_m^{(2)}(x) e^{\lambda _2 x}). y=xk(Rl(1)(x)eλ1x+(Rl(2)(x)eλ2x+Rm(1)(x)eλ1x+(Rm(2)(x)eλ2x).很容易将上式组合成为
      y ∗ = x k e α x [ R n ( 1 ) ( x ) c o s β x + R n ( 2 ) ( x ) s i n β x ] , n = m a x { l , m } , 得 证 . y^* = x^ke^{\alpha x}[R_n^{(1)}(x)cos\beta x + R_n^{(2)}(x)sin\beta x],\quad n = max\{l,m\},得证. y=xkeαx[Rn(1)(x)cosβx+Rn(2)(x)sinβx],n=max{l,m},.
      同样给出简单例子:
      在这里插入图片描述

    3.3.2 拉普拉斯变换法

    \qquad 积分上看Laplace变换是加了条件的Fourier变换,限制了函数的收敛速率 e − s t e^{-st} est。无论是在高阶常微分方程还是常微分方程组的求解过程中都发挥着重要的作用。

    • 定义: L [ f ( x ) ] = ∫ 0 ∞ f ( t ) e − s t d t L[f(x)] = \int_0^{\infty}f(t)e^{-st}dt L[f(x)]=0f(t)estdt
    • 性质1:线性性质, L [ α f ( x ) + β g ( x ) ] = α L [ f ( x ) ] + β L [ g ( x ) ] L[\alpha f(x)+\beta g(x)] =\alpha L[f(x)] +\beta L[g(x)] L[αf(x)+βg(x)]=αL[f(x)]+βL[g(x)]
    • 性质2:原函数的微分性质
      L [ f ′ ( t ) ] = ∫ 0 ∞ f ′ ( t ) e − s t d t L[f'(t)] = \int_0^{\infty}f'(t)e^{-st}dt L[f(t)]=0f(t)estdt = f ( t ) e − s t ∣ 0 ∞ + s ∫ 0 ∞ f ( t ) e − s t d t = f(t)e^{-st}\mid _0^\infty + s\int_0^{\infty}f(t)e^{-st}dt =f(t)est0+s0f(t)estdt = L [ f ( t ) ] − f ( 0 ) = L[f(t)] - f(0) =L[f(t)]f(0)
    • 性质3:像函数的微分性质
      L ′ [ f ( t ) ] = d d s ∫ 0 ∞ f ( t ) e − s t d t L'[f(t)] = \frac{d}{ds}\int_0^{\infty}f(t)e^{-st}dt L[f(t)]=dsd0f(t)estdt = ∫ 0 ∞ − t f ( t ) e − s t d t = \int_0^{\infty}-tf(t)e^{-st}dt =0tf(t)estdt = L [ − t f ( t ) ] =L[-tf(t)] =L[tf(t)]
    • 性质4:位移性质
      已 知 L [ s ] = L [ e a t f ( t ) ] , 已知L[s] = L[e^{at}f(t)], L[s]=L[eatf(t)] 求 L [ e a t f ( t ) ] = ∫ 0 ∞ f ( t ) e a t e − s t d t 求L[e^{at}f(t)] = \int_0^{\infty}f(t)e^{at}e^{-st}dt L[eatf(t)]=0f(t)eatestdt = ∫ 0 ∞ f ( t ) e − ( s − a ) t d t = \int_0^{\infty}f(t)e^{-(s-a)t}dt =0f(t)e(sa)tdt = L [ s − a ] = L[s-a] =L[sa]
      \qquad 利用以上的性质对常微分方程进行拉氏变换,求出y的拉氏变换Y[s],进行相应的拆分,再拉氏逆变换即可得到高阶常微分方程的特解。

    3.4 可降阶的高阶常微分方程

    对于某些可降阶的常微分方程,我们使用变量代换法进行求解。
    能够进行降阶的高阶常微分方程组形会比较简单,仅举例子就不赘述了。

    • 形如 y ( n ) = f ( x ) y^{(n)}= f(x) y(n)=f(x)等式两端直接多次积分。
    • 形如 F ( x , y ( k ) , y ( k + 1 ) , ⋯   , y ( n ) ) = 0 ( 1 ≤ k ≤ n ) F(x,y^{(k)},y^{(k+1)},\cdots,y^{(n)}) = 0(1\leq k\leq n) F(x,y(k),y(k+1),,y(n))=0(1kn),令 q = y ( k ) q = y^{(k)} q=y(k),降阶为 y ( n − k ) y^{(n-k)} y(nk)
    • 形如 F ( y , y ′ , y ′ ′ , ⋯   , y ( n ) ) = 0 F(y,y',y'',\cdots,y^{(n)}) = 0 F(y,y,y,,y(n))=0;
      等式不显含x,令 d y d x = p \frac {dy}{dx} =p dxdy=p;同理可以求出 d 2 y d x 2 = p d p d y , d 3 y d x 3 = p 2 d 2 p d y 2 + p ( d p d y ) 2 , ⋯ ⋯ \frac {d^2y}{dx^2} =p\frac {dp}{dy},\frac {d^3y}{dx^3} =p^2\frac {d^2p}{dy^2} + p(\frac {dp}{dy})^2,\cdots \cdots dx2d2y=pdydp,dx3d3y=p2dy2d2p+p(dydp)2,
    • 形如 F ( x , y ′ , y ′ ′ , ⋯   , y ( n ) ) = d d x Φ ( x , y ′ , y ′ ′ , ⋯   , y ( n − 1 ) ) = 0 F(x,y',y'',\cdots,y^{(n)}) = \frac d{dx}\Phi (x,y',y'',\cdots,y^{(n-1)}) = 0 F(x,y,y,,y(n))=dxdΦ(x,y,y,,y(n1))=0,等式左边为右边的导数,利用右边积分进行降阶。

    值得注意的是,在求解一阶常微分方程组的时候我们会发现,高阶常微分方程实际上可以转换成一阶常微分方程组。

    展开全文
  • matlab解常微分方程

    千次阅读 2021-03-09 20:19:52
    常微分方程ordinary differential equation的缩写,此种表述方式常见于编程,如MATLAB中Simulink求解器solver已能提供了7种微分方程求解方法:ode45(Dormand-Prince),ode23(Bogacki-Shampine),ode113(Adams),ode...
  • Python求解常微分方程

    千次阅读 2021-02-03 16:31:58
    sympy、numpy、scipy、matplotlib是强大的处理数学问题的库,可以执行积分、求解常微分方程、绘图等功能,其开源免费的优势可以与MATLAB媲美。 一阶常微分方程 from sympy import * f = symbols('f', cls=Function...
  • 常微分方程求解器ODE solver

    千次阅读 2021-04-19 23:15:49
    对于常微分方程的初值问题 {dydt=f(t,y),t∈[a,b]y(a)=y0(1) \left\{ \begin{array}{l} \frac{{dy}}{{dt}} = f(t,y),t \in \left[ {a,b} \right]\\ ...定义一个ODEsolver类来使用各种方法求解常微分方程,使用待
  • 常微分方程解法 (四): Matlab 解法

    万次阅读 多人点赞 2019-04-30 10:33:01
    常微分方程解法求解系列博文: 常微分方程解法 (一): 常微分方程的离散化 :差商近似导数、数值积分方法、Taylor 多项式近似 常微分方程解法 (二): 欧拉(Euler)方法 常微分方程解法 (三): 龙格—库塔...
  • 参考:《高等应用数学问题的MATLAB求解(第四版)》 本文涵盖常微分方程的解析法、数值法(Euler、Runge-Kutta法等)
  • 基于MATLAB的高阶常微分方程求解(附完整代码)

    千次阅读 多人点赞 2022-05-02 22:15:14
    单个高阶常微分方程 例题1 二. 高阶常微分方程组 例题2 三. 刚性微分方程 例题3 例题4 四. 隐式微分方程 例题5 一. 单个高阶常微分方程 一个高阶常微分方程的一般形式如下: 输出变量y(t)的各阶导数...
  • MATLAB解常微分方程式的语法是dsolve('equation','condition'),其中equation代表常微分方程式即y'=g(x,y),且须以Dy代表一阶微分项y' D2y代表二阶微分项y'' ,condition则为初始条件。假设有以下三个一阶常微分...
  • 常微分非线性方程式用scipy.integrate.odeint(func, y0, t, args=())
  • 求解一阶常微分方程的数值方法(单步和多步)。 方法包括: 1.欧拉方法2. 亨氏法3. 四阶 Runge Kutta 方法4. Adams-Bashforth 方法5. Adams-Moulton 方法这些方法通常用于求解 IVP,一阶初始值问题 (IVP) 被定义为...
  • 文章目录(1)偏微分方程的类型(二阶)(2)抛物线型1.显式法2.Crank-Nicholson隐式算法 (3)双曲线型(4)椭圆型 (1)偏微分方程的类型(二阶) a∂2u∂x2+b∂2u∂y∂x+c∂2u∂x2+d∂u∂x+e∂u∂y+fu+g=0a\frac{\partial^2u}{\...
  • matlab常微分方程数值求解

    千次阅读 2020-09-01 21:08:13
    本节将介绍常微分方程初值问题的数值求解,主要内容分为三个部分:常微分方程数值求解的概念、求解函数及刚性问题。 一、常微分方程数值求解的一般概念 首先,凡含有参数,未知函数和未知函数导数 (或微分) 的方程,...
  • 常微分方程的数值求解及python实现,包括欧拉公式的理论推导以及相应的python实现代码
  • 文章目录前言1 常微分方程1.1 常微分方程的概念1.2 常微分方程数值求解的一般概念2 常微分方程数值求解函数3 刚性问题结语 前言 1 常微分方程 1.1 常微分方程的概念 1.2 常微分方程数值求解的一般概念 求解常微分...
  • 一份较为详细的神经常微分方程学习笔记,共计194页的。其主要涉及问题定义、ODE的基础理论、不同的求解方法、实验以及应用。欢迎感兴趣的朋友下载学习
  • 神经网络求解二阶常微分方程

    千次阅读 热门讨论 2020-11-20 15:10:41
    基于相关程序源码,我将他的一阶常微分方程求解扩充到二阶常微分方程求解。并且按照此方法可以求解高阶常微分方程。 理论分析 对于任意一个微分方程,我们都可以用这个方程表示出 求解目的就是找出这样的一个方程:...
  • Maple笔记2--常微分方程求解

    千次阅读 2021-02-07 02:24:03
    来源:网络论坛转载(VB资料库)常微分方程求解微分方程求解是数学研究与应用的一个重点和难点. Maple能够显式或隐式地解析地求解许多微分方程求解.在常微分方程求解器dsolve中使用了一些传统的技术例如laplace变换和...
  • MATLAB-6-5常微分方程数值求解

    千次阅读 2020-06-12 15:29:23
    常微分方程数值求解1. 常微分方程数值求解的一般概念2. 常微分方程数值求解函数3. 常微分方程数值求解函数统一命名格式4. 刚性问题 1. 常微分方程数值求解的一般概念 2. 常微分方程数值求解函数 [t,y]=solver...
  • 常微分方程的数值解法

    千次阅读 2021-05-02 16:03:29
    常微分方程的数值解法 常微分方程的分类 初值问题 {dydx=f(x,y)x∈[a,b]y(x0)=y0 \left\{ \begin{aligned} \frac{dy}{dx}=f(x,y)\quad{x}\in[a,b] \\ y(x_0)=y_0\quad \end{aligned} \right . ⎩⎨⎧​dxdy​=f(x,y)...
  • 【实例简介】和Matlab应用有关的,具体介绍常微分方程的使用和解法,原理性介绍,帮助理解。局部截断误差指的是,按()式计算由到这一步的计算值与精确值之差+。为了估计它,由展开得到的精确值是()、()两式相减(注意到...
  • 定义反应活性S=C/(C+D),求ABCD和S的浓度-时间曲线 列微分方程组: 代码: import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def rxt(Z,t): k1=1.0 k2=1.5 r1 = k1 * Z[0...
  • matlab在常微分方程数值解中的应用.docx MATLAB在常微分方程数值解中的应用摘要】许多现实问题都可以通过微分方程的形式进行表示,传统解微分方程的方法【有近似分析解法、表解法和图解法,这些方法需对其进行大量的...
  • python解常微分方程(组)

    千次阅读 多人点赞 2022-01-12 21:15:46
    本文主要讨论了我最近在sympy和odeint上计算常微分方程的一些经验和理解

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 5,026
精华内容 2,010
关键字:

常微分方程求解定义