精华内容
下载资源
问答
  • 偏微分方程特征线

    万次阅读 多人点赞 2018-03-21 23:13:39
    本文以尽可能清晰、简明的方式来介绍了一阶偏微分方程特征线法。个人认为这是偏微分方程理论中较为简单但事实上又容易让人含糊的一部分内容,因此尝试以自己的文字来做一番介绍。当然,更准确来说其实是笔者自己的...

    特别鸣谢原文作者:苏剑林

    本文以尽可能清晰、简明的方式来介绍了一阶偏微分方程的特征线法。个人认为这是偏微分方程理论中较为简单但事实上又容易让人含糊的一部分内容,因此尝试以自己的文字来做一番介绍。当然,更准确来说其实是笔者自己的备忘。

    拟线性情形 ↺

    一般步骤 ↺

    考虑偏微分方程

    α(x,u)xu=β(x,u)(1) (1) α ( x , u ) ⋅ ∂ ∂ x u = β ( x , u )


    其中 α α 是一个 n n 维向量函数,β是一个标量函数, 是向量的点积, uu(x) u ≡ u ( x ) n n 元函数,x是它的自变量。

    特征线法的思路是,设想 x x 是某个参数 s s 的函数,这时候x(s)实际上就是某条高维曲线的参数方程,也就是所谓的特征线,这样 u u 也成为了参数s的函数。并且我们有

    duds=uxdxds(2) (2) d u d s = ∂ u ∂ x ⋅ d x d s


    对比原来的偏微分方程 (1) ( 1 ) ,我们发现可以让

    dxds=α(x,u)(3) (3) d x d s = α ( x , u )


    那么就有

    duds=β(x,u)(4) (4) d u d s = β ( x , u )


    联合 (3) ( 3 ) (4) ( 4 ) ,我们就得到了一个 n+1 n + 1 个方程组成的常微分方程组了。

    <br/>dxds=α(x,u)duds=β(x,u)(5) (5) { d x d s = α ( x , u ) < b r / > d u d s = β ( x , u )


    由于 s s 只是一个额外引入的变量,所以原则上我们可以解得与s无关的结果

    c=f(x,u)(6) (6) c = f ( x , u )


    其中 c c 是一个 n n 维向量,是该常微分方程组的积分常数,f n n 维向量函数。剩下的就是根据初值条件来确定各个积分常数之间的关系了。当然,如果要求出一个通解表达式,那就是

    (7)G(f(x,u))=0


    其中 G G 是任意的n元函数,原则上我们可以从中求解得到 u u 关于x的表达式。

    简单例子 ↺

    上面的步骤比较理论化,在实际处理问题中,可以更加灵活些。下面我们来求解

    ut+xux=u2,u(x,0)=f(x)(8) (8) ∂ u ∂ t + x ∂ u ∂ x = u 2 , u ( x , 0 ) = f ( x )


    我们得到特征线方程

    dt=dxx=duu2(9) (9) d t = d x x = d u u 2


    求得

    x=C1et,u=1C2t(10) (10) x = C 1 e t , u = 1 C 2 − t


    t=0 t = 0 时,有 x=C1,u=1C2=f(C1) x = C 1 , u = 1 C 2 = f ( C 1 ) ,从而解得 C2=1f(C1) C 2 = 1 f ( C 1 ) ,因为我们有 C1=xet,C2=u1+t C 1 = x e − t , C 2 = u − 1 + t ,所以代入得到

    u1+t=1f(xet)(11) (11) u − 1 + t = 1 f ( x e − t )




    u=f(xet)1t×f(xet)(12) (12) u = f ( x e − t ) 1 − t × f ( x e − t )

    一点讨论 ↺

    特征线究竟是什么含义呢?对于初学者来说,上述过程可能像变魔术一样,先求解出常数,然后消去常数,好像把握不到主干。这就是笔者一开始学习特征线法的困惑。

    事实上,我们可以这样认为,特征线已经是偏微分方程的解了,只不过它是在解上的一条线,而完整来说解应该是一个高维的曲面。显然,点动成线、线动成面,想办法让这些线动起来,那么就可以得到这个曲面的方程了,也就是说,各个积分常数 c c 要“动起来”

    当然,它们不能毫无约束地乱动,随意动的话,可能就覆盖了整个空间了,该怎么动是取决于初值条件的,因此我们根据初值条件来确定积分常数之间的约束关系。确定完后,我们已经得到了这个曲面的参数方程,从求解的角度,我们没有必要消去各个常数,但很多时候我们更喜欢显式解,所以我们想办法消去常数。

    整个过程的逻辑大概就是这样了。

    一般情形 ↺

    一般的教材对特征线法的介绍仅限于拟线性偏微分方程,事实上,对于一般的一阶偏微分方程

    F(x,u,ux)=0(13) (13) F ( x , u , ∂ u ∂ x ) = 0


    特征线法也是适用的,其中 F F 是任意多元函数。

    这部分工作主要参考自英文维基百科:
    https://en.wikipedia.org/wiki/Method_of_characteristics

    推导

    为此,我们先记

    (14)p=ux


    然后对 F(x,u,p)=0 F ( x , u , p ) = 0 两边求导,得到

    0=<br/>Fxdxds+Fuuxdxds+Fpdpds=(Fx+Fup)dxds+Fpdpds<br/>(15) (15) 0 = ∂ F ∂ x ⋅ d x d s + ∂ F ∂ u ∂ u ∂ x ⋅ d x d s + ∂ F ∂ p ⋅ d p d s < b r / > = ( ∂ F ∂ x + ∂ F ∂ u p ) ⋅ d x d s + ∂ F ∂ p ⋅ d p d s < b r / >


    可以发现,上面是两组向量的点积,它们加起来为0,那么一个有意思的解决方法就是让

    dxds=Fp,dpds=FxFup(16) (16) d x d s = ∂ F ∂ p , d p d s = − ∂ F ∂ x − ∂ F ∂ u p


    并且我们还有

    duds=uxdxds=pFp(17) (17) d u d s = ∂ u ∂ x ⋅ d x d s = p ⋅ ∂ F ∂ p


    联合起来我们就得到常微分方程组

    <br/><br/><br/>dxds=Fpdpds=FxFupduds=pFpF(x,u,p)=0(18) (18) { d x d s = ∂ F ∂ p < b r / > d p d s = − ∂ F ∂ x − ∂ F ∂ u p < b r / > d u d s = p ⋅ ∂ F ∂ p < b r / > F ( x , u , p ) = 0


    接下来的步骤就跟拟线性情形的没什么差别了,只不过多引入了 n n 个变量p。解出这个方程来,然后得到 2n 2 n 个与 s s 积分常数的积分常数,根据初值条件来找出常数之间的关系。不同的是,因为多了n个变量 p p ,所以要对初值条件也考虑偏导数,这导致了求解过程更加复杂。请看下面的例子。

    又一个例子 ↺

    ut=(ux)2,u(x,0)=f(x)(19) (19) ∂ u ∂ t = ( ∂ u ∂ x ) 2 , u ( x , 0 ) = f ( x )


    也就是 0=F(p)=ptp2x 0 = F ( p ) = p t − p x 2 ,于是根据 (18) ( 18 ) 式我们得到特征线方程

    <br/><br/>dtds=1,dxds=2pxdptds=0,dpxds=0duds=pt2p2x=p2x(20) (20) { d t d s = 1 , d x d s = − 2 p x < b r / > d p t d s = 0 , d p x d s = 0 < b r / > d u d s = p t − 2 p x 2 = − p x 2

    直接取 s=t s = t ,然后可以得到 pt=C1,px=C2 p t = C 1 , p x = C 2 都是常数,从而 x=2C2t+C3,u=C22t+C4 x = − 2 C 2 t + C 3 , u = − C 2 2 t + C 4 。接着根据初值条件,当 t=0 t = 0 时,有 x=C3,u=C4 x = C 3 , u = C 4 ,这表明 C4=f(C3) C 4 = f ( C 3 ) ,代入后有

    u=C22t+f(x+2C2t)(21) (21) u = − C 2 2 t + f ( x + 2 C 2 t )


    然后对初值条件的 x x 变量求偏导,那么当t=0时有

    C2=px=f(x)=f(C3)(22) (22) C 2 = p x = f ′ ( x ) = f ′ ( C 3 )


    注意 f f 是事先给定的初值函数,因此上式实际上是一个代数方程,所以我们有

    (23)u=t×f(C3)2+f(x+2t×f(C3))<br/>=t×f(C3)2+f(C3)<br/>


    并且

    x=2C2t+C3=2t×f(C3)+C3(24) (24) x = − 2 C 2 t + C 3 = − 2 t × f ′ ( C 3 ) + C 3


    从中反解出 C3 C 3 并代入前一式就可以得到完整解。比如当 f(x)=x2 f ( x ) = x 2 时,解得

    C3=x14t(25) (25) C 3 = x 1 − 4 t


    代入得到

    u=x214t(26) (26) u = x 2 1 − 4 t

    自上而下的过程 ↺

    从形式上来看,一般情形的特征线法跟拟线性情形的特征线法差别很大。那么一般情形的特征线法能不能退化到拟线性情形?

    事实上,将 F=αpβ F = α ⋅ p − β 代入 (18) ( 18 ) 式,就可以得到

    <br/>dxds=αduds=pα(27) (27) { d x d s = α < b r / > d u d s = p ⋅ α


    关于 p p 的方程就不必写出来了,因为我们有 pα=β p ⋅ α = β ,代入上式就可以得到封闭的方程组。这就退化为 (5) ( 5 ) 式了。

    另外,一种特别简单的情形是 F F 仅仅是p的函数,这时候 Fx ∂ F ∂ x Fu ∂ F ∂ u 都是0,那么特征线方程中 p p 是常数,从而 x x u u 都是关于s的线性函数,整个方程都是完全可解的。之后,就变成了一个纯粹的代数方程问题了。

    方程组的情形 ↺

    上述特征线技巧能否推广到一阶偏微分方程组?一般是不可行的,因为解决了一阶偏微分方程组就相当于解决了任意阶的偏微分方程了,显然我们没有看到这样的工作(要是可能的话,早就被人做了)。

    不过,要是方程组中偏导算子的部分是一致的话,那么特征线法还是可用的。具体来说,考虑偏微分方程

    (α(x,u)x)u=β(x,u)(28) (28) ( α ( x , u ) ⋅ ∂ ∂ x ) u = β ( x , u )


    这时候 u u 也是一个向量,但是左边的偏导数算子是共享的,右边的 β β 的各个分量可以不一样。这时候也可以用特征线法得到

    <br/>dxds=α(x,u)duds=β(x,u)(29) (29) { d x d s = α ( x , u ) < b r / > d u d s = β ( x , u )

    当然,这只是原来方法的平凡推广。

    原文地址

    展开全文
  • 偏微分方程(PDE) 偏微分方程理论: 物理/工程问题————翻译(建模)/物理工程规律————》数学问题(PDE) 物理/工程问题————求解/数学理论————》数学结果 物理/工程问题————分析————》...

    1. 简介

    微分方程:描述自然界中存在的物理现象和普遍规律。

    • 常微分方程(ODE)
    • 偏微分方程(PDE)

    偏微分方程理论:

    • 物理/工程问题————翻译(建模)/物理工程规律————》数学问题(PDE)
    • 物理/工程问题————求解/数学理论————》数学结果
    • 物理/工程问题————分析————》数学公式/物理意义

    偏微分方程的基本概念:

    • 定义:未知函数及其偏导数所满足的方程;F(x, u(x), Du, D2u,…, Dnu) = 0;
    • 阶数:偏微分方程中偏导数的最高阶数,有n阶就为n。
    • 线性偏微分方程:方程中的未知函数或其偏导数的项都是一次项
    • 齐次方程:每项都含有未知函数或未知函数(定义的u = u(x, y))的偏导数。

    导出微分方程的基本步骤:

    1. 明确研究的物理量
    2. 明确物理量遵守的物理规律
    3. 按物理规律写出微分方程(泛定方程)
    4. 微元法:划出一个微元,分析临近部分和ta的相互作用

    微分方程解决的主要问题:

    • 描述对象特征随时间变化的过程
    • 分析对象特征的变化规律
    • 预报对象特征的未来形态
    • 研究控制对象特征的手段

    常用的物理定律概述

    1. 牛顿第二定律:F = ma;
    2. 胡克定律,在弹性限度内,弹性体的张应力和弹性体的形变量成正比:张应力 = 杨氏模量 * 相对伸长,
    3. 热传导的傅立叶定律:在dt时间内,通过面积元dS流入一个小体积元的热量dQ,与沿着面积元法线方向的温度梯度成正比,也与dS和dt成正比:(k是导热系数,由材料决定)
    4. 牛顿冷却定律:单位时间内从周围介质,传到边界上单位面积的热量,与表面和外界的温度差成正比:(这里u1是外界媒质的温度,h为比例系数)
    5. 扩散定律(斐克定律):单位时间流过某横截面的杂质量dm与该横截面积S和浓度梯度成正比:(式中D为扩散系数,负号是指杂质浓度在减少)

    2. 栗子——香烟过滤嘴的作用

    问题:

    • 过滤嘴的作用和ta的材料与长度有什么关系;
    • 人体吸入的毒品量和哪些因素有关,因素的影响程度怎么样

    模型分析:

    • 分析吸烟时毒物进入人体的模型,建立吸烟过程模型;
    • 假设一个机器人持续吸烟,并且环境不变

    模型假设(假设的变量都在这儿哦):

    • 烟草长 l1,过滤嘴长 l2,总长度 l = l1 + l2,假设毒品均匀分布;
    • 毒物进入空气和沿香烟穿行的数量比:a’ : a, a’ + a = 1;
    • 未点燃烟草和过滤嘴的吸收率(过滤程度)分别为b和ß;
    • 烟雾沿着香烟的穿行速度为常速v,香烟的燃烧速度为u,v >> u;
    • Q为毒物进入人体的总量
    • q(x, t)为毒物的流量,w(x, t)为毒物密度

    模型建立:

    • 建立毒物进入身体的总量:Q = ƒ0T q(l, t)dt, T = l1 / u;

    • 求q(x, 0) = q(x)——流量守恒

      q(x) - q(x + ∆x) = bq(x)∆∂, 0 ≤ x ≤ l1;

      q(x) - q(x + ∆x) = ßq(x)∆∂, l1 ≤ x ≤ l2;

      (∆∂ = ∆x / v)

      同时假设:

      • q(0) = aH0;(香烟毒物被吸入的量)
      • H0 = uw0;(香烟毒物的减少量)

      可得到

      dq/dx = -b/v * q(x), 0 ≤ x ≤ l1;

      dq/dx = -ß/v * q(x), l1 ≤ x ≤ l2;

      &

      q(x) = aH0e-bx/v, 0 ≤ x ≤ l1;

      q(x) = aH0e-bx/ve-ß(x-l1)/v, l1 ≤ x ≤ l2;

    • 目标是求q(l, t),假设t时刻,香烟燃烧到了x = ut;

      将q(x)中的H0拓展为H(t),x随时间会变成x-ut,l1会变成l1 - ut

      可以求得:

      q(x, t) = aH(t)e-b(x-ut)/v,ut ≤ x ≤ l1;

      q(x, t) = aH(t)e-b(l1-ut)/ve-ß(x-l1)/v,l1 ≤ x ≤ l;

    • 求w(ut, t)求∆t内毒物密度的增量。

      w(x, t + ∆t) - w(x, t) = b(q(x, t))/v * ∆t(单位长度烟雾被吸收的部分)

      我们已知:

      • q(x, t) = aH(t)e-b(x-ut)/v;
      • H(t) = uw(ut, t);(H0的拓展公式)

      结果可得到:

      w(ut, t) = w0/a’ * (1 - ae-a’but/v), a’ = 1 - a;

    • 计算Q,吸入一致烟毒物进入人体的总量;

      根据求出来的 w(x, t) 和 q(l, t) 可以代入Q = ƒ0T q(l, t)dt, T = l1 / u;求出现在的Q = aM e-ßl2/v µ®, r = a’bl1/v, µ® = 1 - e-r / r;

    结果分析:

    1. Q和a,M成正比,aM是毒物集中在x = l处的吸入量;

    2. e-ßl2/v是过滤嘴的因素;

    3. µ®是烟草的吸收作用;

    4. 判断与不带过滤嘴的香烟比较

      Q1和Q2的差别为b和ß;

      我们已知:

      • 带滤嘴:Q1 = (aw0v / a’b) e-ßl2/v(1 - e-a’bl1/v);
      • 不带滤嘴:Q~12 = (aw0v / a’b) e-bl2/v(1 - e-a’bl1/v);

      所以明显是滤嘴提高了吸收能力

    特点:

    • 引入两个基本函数:流量q(x,t)和密度w(x,t),运用物理学的守恒定律建立微分方程,构造动态模型。

    3. 偏微方程的导出

    3.1. 波动方程的导出

    • 传输线方程(电报方程 )
    • 均匀薄膜的微小横振动方程
    • 流体力学与声学方程
    • 电磁波方程

    ta们的物理本质根本不同,但这些方程的数学形式与弦振动方程和杆纵振动方程完全一样:

    3.2. 扩散方程的导出

    3.2.1. 细杆热传导

    现象描述:导热细杆各点的温度是不均匀的,热量由高到低传导。

    目的:求出细杆中温度的变化方程。

    定律:

    • 热传导的傅立叶定律:q = -k∆u(单位时间内流过单位时间的热量q与温度的梯度成正比,k为热传导系数);
    • 热量守恒定律:热量变化 = 边界流入 + 热源放出;

    第一种情况(系统无热源):

    (热传导仅由物体内部温度不均引发)

    u(x, t)为x处t时刻的温度。

    一维情况下热传导的傅立叶定律有: q = -k∂u/∂x;(q为热流强度)

    在x方向上(微元法):

    • dt时间流入左表面的热量为:q|x Adt;
    • dt时间流出右表面的热量为:q|x+dx Adt;
    • 所以净流入为:q|x Adt - q|x+dx Adt = -∂q/∂xAdxdt =(然后把q换掉再整合)= kuxxAdxdt;

    因为热量Q = c(øAdx)[u|t - u|t+dt] = c(øAdx)utdt;(假设ø为密度)

    所以更具热量守恒定律得到kuxxAdxdt = c(øAdx)utdt;

    结果:ut - a2uxx = 0, a2 = k/cp(这就得到了均匀物体内部无热源时的热传导方程)

    第二种情况(系统内有热源):

    比第一问不过就是多个热源,也就是加一个热量呗(设F(x, t)为单位时间,单位体积内产生的热量:kuxxAdxdt + F(x, t)Adxdt = c(øAdx)utdt;

    得到:ut - a2uxx = F(x, t) / cp = f(x, t);(ƒ(x, t)为热源强度)

    三维情况下就是:ut - a2∆u = f(x, t);

    3.2.2. 扩散问题

    扩散现象:系统浓度不均用时,物质从高浓度转移向低浓度的现象。

    目的:建立空间各点浓度u(x, y, z, t)的方程。

    物理规律:扩散定律和粒子数守恒定律

    • 扩散定律(裴克定律):单位时间内流过单位面积的分子数或质量,与浓度的梯度成正比:(D为扩散系数)

      沿着n方向的大小:

    • 粒子数守恒定律:单位时间流入一定体积粒子数流出同一体积粒子数的差,等于该体积单位时间内粒子数的增加量。

      同样是微元法,划出一个小立方体v分析浓度变化规律。

      • 计算单位时间沿x-方向的净流入量(负号是表示和浓度梯度的方向相反):(利用两个平面的流入流出的积分)
      • 同理可以求出y方向和z方向的净流入量:
      • 所以总的净流入量为
      • 单位时间的粒子增加数为:
      • 根据粒子数守恒定律就可以得到:

      我们在前面得到过这样的结果:

      所以我们将扩散定律代入就可以得到三维扩散方程

      我们484看到了吗有很多D,如果扩散是均匀的,D就是一个常数,我们可以令D = a2,则有

      那如果这个空间里面有扩散源,也就是扩散是从这儿来的,那么我们可以将ta们相等得到公式:

    3.3. 稳定场方程

    3.3.1. 热传导方程

    如果边界条件与热源内不随时间变化,一定时间后,内部温度会达到稳态。

    u = u(x, y, z), ∂u/∂t = 0;

    so

    ut - a2∆u = f(x, y, z, t) -> ∆u = -f(x, y, z)泊松方程

    3.3.2. 静电场下的泊松方程和Laplace方程

    从高斯定理出发,我们可以结合扩散定律简单的得到:

    我们可以得到

    还因为

    我们可以得到两个方程:

    • 泊松方程:
    • Laplace方程(当的时候):

    4. PDE导出总结

    系统的物理状态除了取决于自己状态,还取决于系统环境的状态。

    • 物理、工程在t>0时刻的系统环境,在数学中称为边界条件;
    • 物理、工程在t=0时刻的系统状态,在数学中称为初始条件;

    定解条件:是边界条件和初始条件的总体,反映了问题的个性。

    泛定方程:不带边界和初始条件的方程,反映问题的个性。

    4.1. 初始条件 ——描述系统的初始状态

    • 波动方程:含有时间的二阶导数,所以需要两个初始条件
      1. u|t=0 = u(x) 系统各点的初始位移
      2. ∂u/∂t|t=0 = v(x) 系统各点的初始速度
    • 热传导方程:含有时间的一阶导数
      1. u(x, t)|t=0 = u(x)初始时刻的温度分布
    • 泊松方程和拉普拉斯方程:不含时间的倒数,不需要初始条件

    注意:初始条件给的是初始状态下物理量的分布,而不是指某一位置

    4.2. 边界条件

    未知函数在边界满足条件:

    1. 第一类Dirchlet边界条件:已知未知函数在边界上的函数值:

      基本形式:u(x, y, z, t)| = ƒ(M, t);

    2. 第二类Neumann:已知未知函数在边界上法线方向的导数值;
      基本形式:∂u(x, t)/∂n|x0 = ƒ(x0, t);

    3. 第三类Robin混合边界条件:混合牛顿冷却定律、傅立叶实验定律(一维);

      牛顿冷却定律:

      q = h(u - ø)n,u为物体表面温度,ø是周围介质温度,h是热交换系数。一维条件下简化为q = h(u - ø)。

      牛顿冷却定律可以得出流出热量与外界温差成正比。


      傅立叶实验定律:

      q|x=a = -k∂u/∂x|x=a

      基本形式:(u + H∂u/∂n)| = ƒ(x0, y0, z0, t);

    4.3. 小结

    5. 偏微分方程的求解

    5.1. 达朗贝尔公式:行波解

    基本思想:

    1. 求PDE的通解;
    2. 用定解求特解

    关键步骤:通过变量变换,将波动方程化为齐次二阶偏微分方程。(所以也叫做波动方程的初始问题,或者柯西问题)

    求解问题:

    我们根据式1可以拆分得到,变换变量也就可以转换为

    通过复合求导我们可以得到:
    ——》

    一些不要的东西删了:

    求通解:

    对两个创建的变量进行积分:

    积分得到

    再积分得到

    我们的通解就是

    其中ƒ1, ƒ2 是连续可微的一元函数。

    ƒ1(x - at)和ƒ2(x + at)的意义?

    • x - at为正方向运动的行波
    • x + at为反方向运动的行波

    确定待定函数形式,求特解:

    我们已经有初始条件:中的后两条。

    我们可以根据之前的条件可以进行两个公式的化简:

    1. image-20210714150847740
    2. image-20210714150921073

    所以可以得到我们求出的达朗贝尔公式:
    image-20210714151011495

    5.2. 分离变量法

    基本思想:将偏微分方程通过分离变量变成一个常微分方程!

    关键步骤:求特征值问题

    适用问题:有界域上的波动方程、热传导方程、稳定场方程等

    求解问题:

    特征值问题:含有未知常数的常微分方程,求非零解的问题;

    特征值:是方程有非零解的常数值;

    特征函数:和特征值相对于的非零解

    5.3. Fourier变换法(傅立叶变换法)

    基本思想:将偏微分方程通过Fourier变换变成一个常微分方程~

    关键步骤:求常微分方程定解问题和ta的解的方法,Fourier逆变换

    适用问题:无界域上的波动方程、热传导方程等

    求解问题:

    Fourier变换法基本步骤:

    1. 对偏微分方程与初始条件中实行傅立叶变换,转化为常微分方程;
    2. 解常微分方程的定解问题,得到相应的傅立叶变换式;
    3. 对该式进行逆傅立叶变换,求的原问题的解
    展开全文
  • Matlab偏微分方程快速上手:使用pdetool工具箱求解二维偏微分方程,适用于数学建模、数学实验,简单的偏微分方程数值计算与工程问题。

    注:本人使用MatlabR2020a版本。

    1.pdetoolbox的调用

    打开MatlabR2020a,在命令行键入pdetool,进入pdetoolbox。
    输入pdetool进入pdetoolbox

    2.绘制定解区域(解的定义域)

    由图形界面可知,解的定义域是 x , y x,y x,y二维坐标构成的平面空间。我们必须设置自己的定解区域,才能定义自己的方程:
    导航栏下方的前5个按钮,分别对应绘制矩形求解区域、绘制按中心生成的矩形求解区域、绘制椭圆形(圆形)定解区域、绘制按中心生成的椭圆形(圆形)定解区域、绘制多边形求解区域。使用时,只需要点击后在绘图区域拖拽(多边形除外,多边形区域是在绘图区域点点以确定顶点),就可以生成定解区域了。
    这是前五个按钮
    上面这是前五个按钮。
    在这里,作者随意绘制了一个椭圆形区域,和一个矩形区域。操作的时候用鼠标拖动操作柄拖拽就可以。(也可以画好几个叠起来)
    在这里,作者随意绘制了一个椭圆形区域,和一个矩形区域
    真的是“随便”画一个就可以哦,因为在matlab下,求解区域的位置坐标精度达到了 1 0 − 16 10^{-16} 1016左右,手动画几乎不可能画准。所以下一步教大家怎么细致地调节边界的坐标。

    3.手动调整定解区域的大小

    双击刚刚绘制好的区域,弹出一个对话框,里面是我们的定解区域的边界坐标信息(注意不全是坐标),我们可以在这里手动调整定解区域的位置(以矩形区域为例):
    刚刚打开时的样子
    这就是刚刚打开时的样子。因为这个区域是作者随便画的,所以坐标信息就像随机数一样。下面我们输入精确的数值:Left: -1, Bottom: -1, Width: 2, Height: 2, Name 就用默认的就好。
    输入参数

    参数的意义:Left:左边界的坐标(x左),Bottom:底边界的坐标(y底),Width:区域宽度,Height:区域高度。这样就得到了 x ∈ [ − 1 , 1 ] , y ∈ [ − 1 , 1 ] x\in[-1,1],y\in[-1,1] x[1,1],y[1,1]的矩形求解区域。调整好的求解区域显示效果如下。绘制好的区域

    4.调整绘图窗口的显示区域(调整显示坐标限)

    有时候我们会发现我们的定解区域太大了,绘图窗口显示不下;或者定解区域太小了,看上去非常不协调。上一个例子中作者的纵坐标显然非常吻合,但是横坐标多出来了(显示了左右两边的白框),那么我强烈建议大家调整完定解区域的坐标以后,再调整一下绘图窗口的显示区域。
    点击导航栏Options,再点击Axes Limits…(意为调整坐标限),可以手动设置坐标限。这里作者勾选了Auto,这样matlab将自动帮我们调整坐标限,使得定解区域位于界面中央。Options还有其他操作,大家可以自己尝试一下,这里就不介绍了。调整坐标限

    调整后的效果如下。

    调整效果

    5.确定边界条件

    点击导航栏下方第6个按钮( ∂ Ω \partial \Omega Ω,意为 Ω \Omega Ω的边界条件),它用来显示边界。下面仍然以矩形区域为例。

    这个按钮
    双击边界

    此时显示了4个边界。双击任意一个边界,会弹出边界条件对话框,可以在这里随意设置边界条件。

    对话框

    在这里可以设置Dirichlet边界条件、Neumann边界条件和Robin边界条件(分别是第一、第二、第三边界条件),其中Robin边界条件和Neumann边界条件集成到一起了。按提示输入对应的系数就可。

    在这里作者使用了如下的边界条件:
    x = ± 1 , u = 0 ; y = ± 1 , ∂ u ∂ n = 0. x=\pm1,u=0;y=\pm1,\frac{\partial u}{\partial n }=0. x=±1,u=0;y=±1,nu=0.
    如果用了Dirichlet边界条件,边界将显示为红色;Neumann、Robin边界条件将显示蓝色,效果如下。

    边界条件

    6.确定偏微分方程的形式

    点击第7个图标(显示PDE字样),按提示输入偏微分方程的系数即可。在这里笔者求解波动方程: ∂ 2 u ∂ 2 t = ∇ u . \frac{\partial^2 u}{\partial^2 t}=\nabla u. 2t2u=u.

    在这里插入图片描述

    工具箱提供的方程通式如下:
    1.椭圆型Elliptic,通用数学形式为 − ∇ ⋅ ( c ∇ u ) + a u = f ; -\nabla \cdot(c\nabla u)+au=f ; (cu)+au=f;
    2.抛物型Parabolic,通用数学形式为 d ∂ u ∂ t − ∇ ⋅ ( c ∇ u ) + a u = f ; d\frac{\partial u}{\partial t}-\nabla \cdot(c\nabla u)+au=f ; dtu(cu)+au=f;
    3.双曲型Hyperbolic,通用数学形式为 d ∂ 2 u ∂ 2 t − ∇ ⋅ ( c ∇ u ) + a u = f ; d\frac{\partial^2 u}{\partial^2 t}-\nabla \cdot(c\nabla u)+au=f ; d2t2u(cu)+au=f;
    4.特征值方程Eigenmodes,若 λ \lambda λ为特征值,则数学形式为 − ∇ ⋅ ( c ∇ u ) + a u = λ d u . -\nabla \cdot(c\nabla u)+au=\lambda d u . (cu)+au=λdu.

    也可以自行指定求解的方程类型,比如比较常见的热传导、扩散等方程,可以在下面图示的下拉菜单中选择,但是仍然要按上面讲的方法手动设置系数。

    手动设置方程类型

    7.三角剖分

    由于Matlab pdetoolbox使用有限元方法求解,所以需要三角剖分。点击第8个图标(1个三角形图样)可以初始化剖分,点击第9个图标(4个三角形图样)可以增加剖分密度,这样可以提高计算精度,但是密度过高内存可能会爆掉,使用要谨慎。
    三角剖分

    8.设置初始条件,准备求解

    点击导航栏Solve,再点击Parameters…,进入求解参数设置器。在这里,第一行Time我们可以设置 t t t的求解范围及步长,默认情况下是不显示步长的(默认显示0:10意为从0求解到10,步长为1),我们按照Maltab等差数列的生成方法 a:j:b 就可以设置时间步长j了。

    第二行u(t0)、第三行u’(t0)表示 t 0 t_0 t0时刻的两个初始条件。这里作者使用了如下的初始条件。

    初始条件
    第四行和第五行表示相对容差和绝对容差,笔者查看了Matlab帮助中心,大概了解到这两个参数似乎与浮点数0的截断精度有关,太小的话会延长计算时间,如果你想了解更多,笔者把链接提供上来Absolute tolerance - MATLAB & Simulink - MathWorks 中国,假如我们对计算精度没有要求的话,使用默认值就可以了。这里笔者为了演示使用了0.001和0.0001。如果想跟着一起做,那么笔者把方程的代码也放上来:第一个是atan(cos(pi/2*x)),第二个是3*sin(pi*x).*exp(cos(pi*y))

    9.求解

    点击导航栏下方按钮(一个“ = = =”字样的按钮,就是增加三角剖分密度右边那个按钮),这个按钮表示开始求解。如果求解完成的话会显示这个图。在这里可以点击“放大镜”按钮寻找感兴趣的区域放大来观察细节(放大之后想要缩小就要用上面步骤4的方法重新设置坐标限了,没有找到缩小的快捷键)。求解结果
    它直接显示了 t = 10 t=10 t=10 u u u 求 解 区 域 Ω 求解区域\Omega Ω的图像。这样的输出缺乏直观性,我们点击导航栏下方一个长得像matlab的logo的按钮(就是“ = = =”按钮右边那个),调整绘图格式。

    输出格式窗口
    这个窗口有许多功能,作者就不再一一详述了。大家可以自行调试。比较常用的有“Contour”绘制等高线图,“Arrows”绘制向量场,“Height(3-D plot)”按3D模式输出(这个比较常用),“Animation”按动画形式输出(2D\3D都支持),此选项勾选后右边的Option选项会变亮,我们可以点进去在里面设置1秒显示的帧数、重复播放的次数。这里作者按照25阶变色、‘jet’ Colormap、向量场为 − ∇ u -\nabla u u的形式静态输出 t = 0.5 t=0.5 t=0.5时的结果,如下图所示。
    在这里插入图片描述

    这就是matlab pdetool工具箱的主要使用方法,本人也是小白一枚,所以欢迎大家批评指正,可以在评论区留下你的想法。

    参考:偏微分方程(姜礼尚《数学物理方程讲义》第三版)(更新完毕,附课件)来自于西北大学数统学院的马老师的数理方程视频课,是按数学系的讲法讲的数理方程,里面有那么两三个视频是讲如何用Matlab求解偏微分方程,如果你懂偏微分方程的话进去听一遍就会了。
    感觉应该是疫情期间这位老师的网课视频?那么西北大学的学生们也太幸福了,因为马老师讲的真的很好!人也很好玩哈哈,还去b站里别的老师的微分几何课程下面评论,正巧那个被评论的老师的同学就在b站讲拓扑哈哈,那个老师讲的也特别细我还给听完了(浙江理工庄老师)。扯远了!但是还是安利!!!

    展开全文
  • 高等数学 偏微分方程 波动方程


    偏微分方程(Partial Differential Equation I)
    偏微分方程(Partial Differential Equation II)
    偏微分方程(Partial Differential Equation III)
    偏微分方程(Partial Differential Equation IV)


    参考文献:

    《数学物理方程》| 季孝达
    《数学物理方法》| 吴崇试
    《数学物理方法》| 梁昆淼
    MOOC北京大学《数学物理方法》| 吴崇试 、高春媛

    偏微分方程的定解问题

    偏微分方程 (Partial Differential Equations, PDE) ,是指含有多元未知函数及其若干阶偏导数的微分方程。
    F ( u , x 1 , x 2 , ⋯   , ∂ u ∂ x 1 , ∂ u ∂ x 2 , ⋯   , ∂ 2 u ∂ x 1 2 , ⋯   ) = 0 (1.1) F(u,x_1,x_2,\cdots,\cfrac{∂u}{∂x_1},\cfrac{∂u}{∂x_2},\cdots,\cfrac{∂^2u}{∂x_1^2},\cdots)=0\tag{1.1} F(u,x1,x2,,x1u,x2u,,x122u,)=0(1.1)
    其中最高阶偏导数的阶称为偏微分方程的阶。

    在此引入一些重要的微分算子

    Hamilton 算子: ∇ = ( ∂ ∂ x , ∂ ∂ y , ∂ ∂ z ) ∇=(\cfrac{∂}{∂x},\cfrac{∂}{∂y},\cfrac{∂}{∂z}) =(x,y,z)
    Laplace算子 Δ = ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 = ∇ ⋅ ∇ = ∇ 2 Δ=\cfrac{∂^2}{∂x^2}+\cfrac{∂^2}{∂y^2}+\cfrac{∂^2}{∂z^2}=∇\cdot∇=∇^2 Δ=x22+y22+z22==2

    定解问题及适定性

    数学物理方程 (Equation of Mathematical Physics) ,通常指从物理学及其他各门自然科学、技术科学中产生的偏微分方程。数学物理方程有时还包括常微分方程和积分方程。本课程将着重讨论三类基本的二阶线性偏微分方程。

    (1) 弦的横振动方程:未知函数 u ( x , t ) u(x,t) u(x,t) 表示坐标 x x x 处弦的横向位移
    u t t − a 2 u x x = f ( x , t ) u_{tt}-a^2u_{xx}=f(x,t) utta2uxx=f(x,t)
    其中 a = T ρ , f = F ρ a=\sqrt{\cfrac{T}{\rho}},f=\cfrac{F}{\rho} a=ρT ,f=ρF ρ \rho ρ 为弦的线密度, T T T 为弦的切应力, F F F 表示单位长度受到的横向力。

    杆的纵振动方程:未知函数 u ( x , t ) u(x,t) u(x,t) 表示截面相对于平衡位置 x x x 处的纵向位移
    u t t − a 2 u x x = f ( x , t ) u_{tt}-a^2u_{xx}=f(x,t) utta2uxx=f(x,t)
    其中 a = E ρ , f = F ρ a=\sqrt{\cfrac{E}{\rho}},f=\cfrac{F}{\rho} a=ρE ,f=ρF ρ \rho ρ 为杆的线密度, E E E 为杆的Young模量, F F F 表示单位长度受到的纵向力

    更一般的,三维空间的波动方程 (wave equation)是
    u t t − a 2 Δ u = f ( r , t ) (1.2) u_{tt}-a^2Δu=f(\mathbf r,t)\tag{1.2} utta2Δu=f(r,t)(1.2)

    (2) 热传导方程:(heat equation)未知函数 u ( r , t ) u(\mathbf r,t) u(r,t) 表示温度
    u t − a 2 Δ u = f ( r , t ) (1.3) u_t-a^2 Δu=f(\mathbf r,t) \tag{1.3} uta2Δu=f(r,t)(1.3)
    其中 κ ρ c , F ρ c \cfrac{\kappa}{\rho c},\cfrac{F}{\rho c} ρcκ,ρcF κ \kappa κ 称为热传导系数, ρ \rho ρ 是介质的密度, c c c 是比热容,与介质的质料有关。 F ( x , y , z , t ) F(x,y,z,t) F(x,y,z,t) 表示介质中有热源,单位时间内在单位体积介质中产生的热量

    (3) 如果热传导方程热源不随时间变化,当温度达到稳定时 ( u t = 0 ) (u_t=0) (ut=0) ,温度分布满足 Poisson方程
    Δ u = − f / a 2 (1.4) Δu=-f/a^2\tag{1.4} Δu=f/a2(1.4)
    如果没有热源 ( f = 0 ) (f=0) (f=0) ,则有 Laplace方程(也称调和方程)
    Δ u = 0 (1.5) Δu=0\tag{1.5} Δu=0(1.5)
    静电场的电势 u ( r ) u(\mathbf r) u(r)也满足 Poisson方程
    Δ u = − ρ / ϵ 0 Δu=-\rho/\epsilon_0 Δu=ρ/ϵ0
    其中 ρ \rho ρ 为电荷密度, ϵ 0 \epsilon_0 ϵ0 为真空介电常数。

    (4) 如果波动方程 u ( r , t ) u(\mathbf r,t) u(r,t) 随时间周期性变化,频率为 ω \omega ω,则 u ( r , t ) = v ( r ) e − i ω t u(\mathbf r,t)=v(\mathbf r)e^{-i\omega t} u(r,t)=v(r)eiωt,则 v ( r ) v(\mathbf r) v(r) 满足 Helmholtz 方程
    Δ v + k 2 v = 0 (1.6) Δv+k^2v=0\tag{1.6} Δv+k2v=0(1.6)
    其中 k = ω / a k=\omega/a k=ω/a 称为波数。

    通解和特解:如果多元函数 u u u 具有偏微分方程中出现的各阶连续偏导数,并使方程 恒成立,则称此函数为方程的解,也称古典解。和常微分方程类似,称 m m m 阶偏微分方程的含有 m m m 个任意函数的解为通解。通解中的任意函数一旦确定便成为特解。

    示例 1:求解偏微分方程
    u x y = 0 u_{xy}=0 uxy=0
    解:先关于 y y y 积分,得 u x = f ( x ) u_x=f(x) ux=f(x),再关于 x x x 积分,就得到通解
    u = ∫ f ( ξ ) d ξ + f 2 ( η ) = f 1 ( ξ ) + f 2 ( η ) u =\int f(ξ)\mathrm dξ+f_2(η)=f_1(ξ)+f_2(η) u=f(ξ)dξ+f2(η)=f1(ξ)+f2(η)
    其中 f 1 , f 2 f_1,f_2 f1,f2 是任意函数。

    示例 2:二维Laplace方程的通解
    u x x + u y y = 0 u_{xx}+u_{yy}=0 uxx+uyy=0
    引入变换
    { ξ = x + i y η = x − i y \begin{cases} ξ=x+iy \\ η=x-iy \end{cases} {ξ=x+iyη=xiy
    根据复合函数偏导法则,进一步求得
    { u x x = u ξ ξ + 2 u ξ η + u η η u y y = − ( u ξ ξ − 2 u ξ η + u η η ) \begin{cases} u_{xx}=u_{ξξ}+2u_{ξη}+u_{ηη} \\ u_{yy}=-(u_{ξξ}-2u_{ξη}+u_{ηη}) \end{cases} {uxx=uξξ+2uξη+uηηuyy=(uξξ2uξη+uηη)
    原方程变为
    u ξ η = 0 u_{ξη}=0 uξη=0
    于是可求得通解
    u ( x , y ) = f ( x + i y ) + g ( x − i y ) u(x,y)=f(x+iy)+g(x-iy) u(x,y)=f(x+iy)+g(xiy)

    定解条件:通常把反应系统内部作用导出的偏微分方程称为泛定方程。为了完全描写一个实际物理问题的确定解,我们需要在一定的制约条件,称为定解条件。偏微分方程和相应的定解条件就构成一个定解问题。常见的定解条件有以下几类:

    • 初始条件:应该完全描写初始时刻( t = 0 t=0 t=0)介质内部及边界上任意一点的状况。一般的讲,关于时间 t t t m m m 阶偏微分方程,要给出 m m m 个初始条件才能确定一个特解。
      对于波动方程来说,就是初始时刻的位移和速度
      u ∣ t = 0 = ϕ ( r ) , ∂ u ∂ t ∣ t = 0 = ψ ( r ) u|_{t=0}=\phi(\mathbf r),\quad \cfrac{∂u}{∂t}|_{t=0}=\psi(\mathbf r) ut=0=ϕ(r),tut=0=ψ(r)
      对于热传导方程,只用给出初始时刻物体温度的分布状态
      u ∣ t = 0 = ϕ ( r ) u|_{t=0}=\phi(\mathbf r) ut=0=ϕ(r)
      Laplace 方程和Possion 方程都是描述稳恒状态的,与时间无关,所以不需要初始条件。

    • 边界条件:边界条件形式比较多样化,要由具体问题中描述的具体状态决定。总的原则是:边界条件应该完全描写边界上各点在任意时刻的状况。
      对于弦的横振动问题,约束条件通常有以下三种

    1. 如果弦的两端固定,那么边界条件就是
      u ∣ x = 0 = 0 , u ∣ x = l = 0 u|_{x=0}=0,\quad u|_{x=l}=0 ux=0=0,ux=l=0

    2. 如果一端 ( x = 0 ) (x=0) (x=0) 固定,另一端 ( x = l ) (x=l) (x=l) 受位移方向的外力 F ( t ) F(t) F(t) ,可以推导出
      u ∣ x = 0 = 0 , ∂ u ∂ x ∣ x = l = 1 T F ( t ) u|_{x=0}=0,\quad \cfrac{∂u}{∂x}|_{x=l}=\cfrac{1}{T}F(t) ux=0=0,xux=l=T1F(t)
      如果一端固定,另一端外力为 0,即另一端自由,则
      u ∣ x = 0 = 0 , ∂ u ∂ x ∣ x = l = 0 u|_{x=0}=0,\quad \cfrac{∂u}{∂x}|_{x=l}=0 ux=0=0,xux=l=0

    3. 如果外力是由弹簧提供的弹性力,即 F ( t ) = − k u ( l , t ) F(t)=-ku(l,t) F(t)=ku(l,t),其中 k k k 为弹性系数,则
      ( ∂ u ∂ x + σ u ) ∣ x = l = 0 , σ = k T (\cfrac{∂u}{∂x}+\sigma u)|_{x=l}=0,\quad \sigma=\cfrac{k}{T} (xu+σu)x=l=0,σ=Tk

      对于热传导问题,也有类似的边界条件(以 ∂ Ω ∂\Omega Ω 表示区域 Ω \Omega Ω 的边界)

    4. 如果边界温度分布已知,则
      u ∣ r ∈ ∂ Ω = f ( r , t ) u|_{\mathbf r\in∂\Omega}=f(\mathbf r,t) urΩ=f(r,t)

    5. 如果物体边界和周围介质绝热,则
      ∂ u ∂ n ∣ r ∈ ∂ Ω = 0 \cfrac{∂u}{∂n}|_{\mathbf r\in∂\Omega}=0 nurΩ=0

    6. 设周围介质的温度为 u 1 ( r , t ) u_1(\mathbf r,t) u1(r,t) ,物体通过边界与周围有热量交换,则
      ( ∂ u ∂ n + σ u ) ∣ r ∈ ∂ Ω = σ u 1 (\cfrac{∂u}{∂n}+\sigma u)|_{\mathbf r\in∂\Omega}=\sigma u_1 (nu+σu)rΩ=σu1
      其中 σ = h / k , h \sigma=h/k,h σ=h/k,h表示两种物质间的热交换系数。

      常见的边界条件数学上分为三类(以 ∂ Ω ∂\Omega Ω 表示区域 Ω \Omega Ω 的边界, f f f 为已知函数)

    7. 第一类边界条件 (Dircichlet条件):直接给出边界上各点未知函数 u u u 的值
      u ∣ ∂ Ω = f ( r , t ) u|_{∂\Omega}=f(\mathbf r,t) uΩ=f(r,t)

    8. 第二类边界条件 (Neu-mann条件):给出边界外法线方向上方向导数的数值
      ∂ u ∂ n ∣ ∂ Ω = f ( r , t ) \cfrac{∂u}{∂n}|_{∂\Omega}=f(\mathbf r,t) nuΩ=f(r,t)

    9. 第三类边界条件 (Robin条件):给出边界上各点的函数值与外法线方向上方向导数的某一线性组合值
      ( ∂ u ∂ n + σ u ) ∣ ∂ Ω = f ( r , t ) (\cfrac{∂u}{∂n}+\sigma u)|_{∂\Omega}=f(\mathbf r,t) (nu+σu)Ω=f(r,t)

    • 连接条件:从微分方程的推导知道,方程只在空间区域的内部成立,如果在区域内部出现结构上的跃变,那么在这些跃变点(线,面)上微分方程不再成立,因此需要补充上相应的条件,通常称为连接条件
      例如,由两段不同材质组成的弦在 x 0 x_0 x0 处连接,波动方程需分段讨论,且在连接处位移和张力相等,连接条件为 u 1 ( x 0 , t ) = u 2 ( x 0 , t ) , ∂ u 1 ∂ x ∣ x = x 0 − 0 = ∂ u 2 ∂ x ∣ x = x 0 + 0 u_1(x_0,t)=u_2(x_0,t),\cfrac{∂u_1}{∂x}|_{x=x_0-0}=\cfrac{∂u_2}{∂x}|_{x=x_0+0} u1(x0,t)=u2(x0,t),xu1x=x00=xu2x=x0+0

    定解问题的适定性:定解问题来自于实际,它的解也应该回到实际中去。如果一个定解问题的解存在、唯一且稳定,则称该定解问题是适定的 (well-posed) 。稳定性指的是,如果定解条件的数值有细微的改变,解的数值也作细微的改变。

    线性叠加原理

    线性叠加原理:考虑 n n n 个自变量(包括时间 t t t x 1 , x 2 , ⋯   , x n x_1,x_2,\cdots,x_n x1,x2,,xn 的二阶线性偏微分方程
    ∑ j = 1 n ∑ i = 1 n a i j u x i x j + ∑ i = 1 n b i u x i + c u = f \displaystyle\sum_{j=1}^{n}\sum_{i=1}^n a_{ij}u_{x_ix_j}+\sum_{i=1}^{n}b_iu_{x_i}+cu=f j=1ni=1naijuxixj+i=1nbiuxi+cu=f
    其中 a i j , b i , c , f a_{ij},b_i,c,f aij,bi,c,f x 1 , x 2 , ⋯   , x n x_1,x_2,\cdots,x_n x1,x2,,xn 的已知函数, f f f 称为方程的非齐次项。

    引入微分算子
    L = ∑ j = 1 n ∑ i = 1 n a i j ∂ ∂ x i ∂ x j + ∑ i = 1 n b i ∂ ∂ x i + c L=\sum_{j=1}^{n}\sum_{i=1}^n a_{ij}\cfrac{∂}{∂x_i∂x_j} +\sum_{i=1}^{n}b_i\cfrac{∂}{∂x_i}+c L=j=1ni=1naijxixj+i=1nbixi+c
    则可简单的表示为
    L [ u ] = f L[u]=f L[u]=f
    一般的,对任意常数 c 1 , c 2 c_1,c_2 c1,c2 和任意函数 u 1 , u 2 u_1,u_2 u1,u2 ,微分算子满足性质
    L [ c 1 u 1 + c 2 u 2 ] = c 1 L [ u 1 ] + c 2 L [ u 2 ] (2.1) L[c_1u_1+c_2u_2]=c_1L[u_1]+c_2L[u_2]\tag{2.1} L[c1u1+c2u2]=c1L[u1]+c2L[u2](2.1)
    称为线性微分算子。显然, L L L 为线性微分算子。
    对于一般的线性边界条件(包括三类边界条件)也可以写成算子的形式
    L 0 [ u ] = ( α ∂ u ∂ n + β u ) ∣ ∂ Ω = ϕ L_0[u]=(\alpha\cfrac{∂u}{∂n}+\beta u)|_{∂\Omega}=\phi L0[u]=(αnu+βu)Ω=ϕ
    其中 α , β , ϕ \alpha,\beta,\phi α,β,ϕ 是已知函数,容易证明 L 0 L_0 L0 也是线性微分算子。

    叠加原理
    (1) 有限叠加原理:若 u i u_i ui 分别满足方程 L [ u i ] = f i ( i = 1 , 2 , ⋯   , m ) L[u_i]=f_i\quad(i=1,2,\cdots,m) L[ui]=fi(i=1,2,,m) ,则他们的线性组合 u = ∑ i = 1 m c i u i u=\displaystyle\sum_{i=1}^mc_iu_i u=i=1mciui 满足方程
    L [ u ] = ∑ i = 1 m c i f i (2.2) L[u]=\displaystyle\sum_{i=1}^mc_if_i\tag{2.2} L[u]=i=1mcifi(2.2)
    (2) 级数叠加原理:若 u i u_i ui 分别满足方程 L [ u i ] = f i ( i = 1 , 2 , ⋯   ) L[u_i]=f_i\quad(i=1,2,\cdots) L[ui]=fi(i=1,2,) ,则他们的级数 u = ∑ i = 1 ∞ c i u i u=\displaystyle\sum_{i=1}^\infty c_iu_i u=i=1ciui 满足方程
    L [ u ] = ∑ i = 1 ∞ c i f i (2.3) L[u]=\displaystyle\sum_{i=1}^\infty c_if_i\tag{2.3} L[u]=i=1cifi(2.3)
    (3) 积分叠加原理:若 u ( r ; r 0 ) u(\mathbf{r;r_0}) u(r;r0) 满足方程 L [ u ] = f ( r ; r 0 ) L[u]=f(\mathbf{r;r_0}) L[u]=f(r;r0) ,则积分 U ( r ) = ∫ V c ( r 0 ) u ( r ; r 0 ) d r 0 U(\mathbf r)=\displaystyle\int_V c(\mathbf r_0)u(\mathbf{r;r_0})\mathrm d\mathbf r_0 U(r)=Vc(r0)u(r;r0)dr0 满足方程
    L [ U ] = ∫ V c ( r 0 ) f ( r ; r 0 ) d r 0 (2.4) L[U]=\int_Vc(\mathbf r_0)f(\mathbf{r;r_0})\mathrm d\mathbf r_0\tag{2.4} L[U]=Vc(r0)f(r;r0)dr0(2.4)

    一阶(拟)线性偏微分方程

    引例:求解右行单波方程
    u t + a u x = 0 u_{t}+au_{x}=0 ut+aux=0
    的通解,引入变换
    { ξ = x − a t η = x \begin{cases} ξ=x-at \\ η=x \end{cases} {ξ=xatη=x
    根据复合函数偏导法则,进一步求得
    { u t = − a u ξ u x = u ξ + u η \begin{cases} u_{t}=-au_{ξ}\\ u_{x}=u_{ξ}+u_{η} \end{cases} {ut=auξux=uξ+uη
    原方程变为
    u η = 0 u_{η}=0 uη=0
    于是可求得通解
    u ( x , t ) = g ( x − a t ) u(x,t)=g(x-at) u(x,t)=g(xat)

    一阶线性偏微分方程

    n n n 个自变量的一阶线性偏微分方程的一般形式为
    ∑ i = 1 n b i ∂ u ∂ x i + c u = f \sum_{i=1}^{n}b_i\cfrac{∂u}{∂x_i}+cu=f i=1nbixiu+cu=f
    其中 b i , c , f b_i,c,f bi,c,f 是自变量 x 1 , x 2 , ⋯   , x n x_1,x_2,\cdots,x_n x1,x2,,xn 的函数。

    先研究两自变量 x , y x,y x,y 的一阶线性偏微分方程
    a u x + b u y + c u + f = 0 (2.1) au_x+bu_y+cu+f=0\tag{2.1} aux+buy+cu+f=0(2.1)
    其中 a , b , c , f a,b,c,f a,b,c,f 都是自变量 x , y x,y x,y 的函数。
    (1) 若 a ≡ 0 , b ≠ 0 a\equiv0,b\neq0 a0,b=0 方程改写为
    u y + c b u + f b = 0 u_y+\cfrac{c}{b}u+\cfrac{f}{b}=0 uy+bcu+bf=0
    利用一阶线性常微分方程的求解方法可求得通解
    u = e − p ( x , y ) [ ∫ e p ( x , y ) f ( x , y ) b ( x , y ) d y + g ( x ) ] u=e^{-p(x,y)}[\int e^{p(x,y)}\cfrac{f(x,y)}{b(x,y)} \mathrm dy+g(x)] u=ep(x,y)[ep(x,y)b(x,y)f(x,y)dy+g(x)]
    其中 p ( x , y ) = ∫ c ( x , y ) b ( x , y ) d y , g ( x ) p(x,y)=\displaystyle\int\cfrac{c(x,y)}{b(x,y)} \mathrm dy,g(x) p(x,y)=b(x,y)c(x,y)dy,g(x) 是任意函数。
    (2) 若 a ( x , y ) b ( x , y ) ≠ 0 a(x,y)b(x,y)\neq0 a(x,y)b(x,y)=0 不能直接积分求解。我们希望通过适当的自变量变换和未知函数的变换,使方程简化,并在此基础上对方程求解。
    首先作自变量非奇异变换(可逆)
    { ξ = ξ ( x , y ) η = η ( x , y ) \begin{cases} ξ=ξ(x,y) \\ η=η(x,y) \end{cases} {ξ=ξ(x,y)η=η(x,y)
    要求雅可比行列式
    ∣ J ( ξ , η ) ∣ = ∣ ∂ ( ξ , η ) ∂ ( x , y ) ∣ = ∣ ξ x ξ y η x η y ∣ ≠ 0 |J(ξ,η)|=|\cfrac{∂(ξ,η)}{∂(x,y)}|=\begin{vmatrix}ξ_x&ξ_y \\ η_x&η_y\end{vmatrix}\neq0 J(ξ,η)=(x,y)(ξ,η)=ξxηxξyηy=0
    以保证新变量 ξ , η ξ,η ξ,η 相互独立,利用链式法则
    u x = ∂ u ∂ ξ ∂ ξ ∂ x + ∂ u ∂ η ∂ η ∂ x , u y = ∂ u ∂ ξ ∂ ξ ∂ y + ∂ u ∂ η ∂ η ∂ y u_x=\cfrac{∂u}{∂ξ}\cfrac{∂ξ}{∂x}+\cfrac{∂u}{∂η}\cfrac{∂η}{∂x},\quad u_y=\cfrac{∂u}{∂ξ}\cfrac{∂ξ}{∂y}+\cfrac{∂u}{∂η}\cfrac{∂η}{∂y} ux=ξuxξ+ηuxη,uy=ξuyξ+ηuyη
    可得到新的方程
    A u ξ + B u η + C u + F = 0 (3.2) Au_ξ+Bu_η+Cu+F=0\tag{3.2} Auξ+Buη+Cu+F=0(3.2)
    其中
    { A = a ξ x + b ξ y B = a η x + b η y C = c , F = f (3.3) \begin{cases} A=aξ_x+bξ_y \\ B=aη_x+bη_y \\ C=c,F=f \end{cases} \tag{3.3} A=aξx+bξyB=aηx+bηyC=c,F=f(3.3)
    新方程仍然是线性的。从上式看出,如果取 ξ = ξ ( x , y ) ξ=ξ(x,y) ξ=ξ(x,y) 是一阶齐次线性偏微分方程
    a ∂ z ∂ x + b ∂ z ∂ y = 0 (3.4) a\cfrac{∂z}{∂x}+b\cfrac{∂z}{∂y}=0\tag{3.4} axz+byz=0(3.4)
    的解,则此时 A = 0 A=0 A=0,这样方程 (3.2) 得以化简为
    B u η + C u + F = 0 Bu_η+Cu+F=0 Buη+Cu+F=0
    从而可以像第一种类型一样积分求解。
    (3) 偏微分方程 (3.4) 可以做如下几何解释,方程可以改写为向量形式
    ( a , b ) ⋅ ( ∂ z ∂ x , ∂ z ∂ y ) = 0 (3.5) (a,b)\cdot(\cfrac{∂z}{∂x},\cfrac{∂z}{∂y})=0\tag{3.5} (a,b)(xz,yz)=0(3.5)
    方程的解 z = ϕ ( x , y ) z=\phi(x,y) z=ϕ(x,y) 表示空间 x y z xyz xyz 中一张曲面 S : z = ϕ ( x , y ) S:z=\phi(x,y) S:z=ϕ(x,y) ,任取平面 z = c z=c z=c 截得的曲线方程为 L : ϕ ( x , y ) = c L:\phi(x,y)=c L:ϕ(x,y)=c ,如图

    曲线 L L L 上任意一点 P ( x , y ) P(x,y) P(x,y) 的法线方向为 ( ∂ ϕ ∂ x , ∂ ϕ ∂ y ) (\cfrac{∂\phi}{∂x},\cfrac{∂\phi}{∂y}) (xϕ,yϕ) ,式 (3.5) 表明在 P P P 点的向量 ( a , b ) (a,b) (a,b) 与法线方向垂直,即切线方向。 P P P 点切向量微元可以表示为 ( d x , d y ) (\mathrm dx,\mathrm dy) (dx,dy) ,所以存在对应关系
    d x a ( x , y ) = d y b ( x , y ) (3.6) \cfrac{\mathrm dx}{a(x,y)}=\cfrac{\mathrm dy}{b(x,y)}\tag{3.6} a(x,y)dx=b(x,y)dy(3.6)
    由于平面 z = c z=c z=c 的任意性,所截得的 L L L 形成的曲线簇覆盖整个曲面 S S S。所以,只要求出常微分方程 (3.6) 的积分曲线簇 ϕ ( x , y ) = c \phi(x,y)=c ϕ(x,y)=c ,便可知道偏微分方程 (3.4) 对应的解 z = ϕ ( x , y ) z=\phi(x,y) z=ϕ(x,y) 。反之,偏微分方程 (3.4) 的解也是常微分方程 (3.6) 的积分曲线簇。
    证明:假设常微分方程 (3.6) 隐式通解为 ϕ ( x , y ) = c \phi(x,y)=c ϕ(x,y)=c ,则
    d ϕ ( x , y ) = ∂ ϕ ∂ x d x + ∂ ϕ ∂ y d y = 0 \mathrm d\phi(x,y)=\cfrac{∂\phi}{∂x}\mathrm dx+\cfrac{∂\phi}{∂y}\mathrm dy=0 dϕ(x,y)=xϕdx+yϕdy=0
    将常微分方程改写成 d y d x = b ( x , y ) a ( x , y ) \cfrac{\mathrm dy}{\mathrm dx}=\cfrac{b(x,y)}{a(x,y)} dxdy=a(x,y)b(x,y) 带入上式可得
    a ∂ ϕ ∂ x + b ∂ ϕ ∂ y = 0 a\cfrac{∂\phi}{∂x}+b\cfrac{∂\phi}{∂y}=0 axϕ+byϕ=0
    可知 z = ϕ ( x , y ) z=\phi(x,y) z=ϕ(x,y) 是偏微分方程 (3.4) 的一个解。反向推导亦成立。

    (4) 常微分方程 (3.6) 叫做偏微分方程 (3.1) 的特征方程,特征方程的积分曲线叫做特征线。如果求出了积分曲线簇 ϕ ( x , y ) = c \phi(x,y)=c ϕ(x,y)=c ,再任取常数 c c c 使得 ∣ J ( ξ , η ) ∣ ≠ 0 |J(ξ,η)|\neq0 J(ξ,η)=0 ,以此 ξ , η ξ,η ξ,η 作为变量代换,则一阶线性偏微分方程便可求得通解。

    一阶拟线性偏微分方程

    n n n 个自变量的一阶拟线性偏微分方程的一般形式为
    ∑ i = 1 n b i ( x 1 , x 2 , ⋯   , x n , u ) ∂ u ∂ x i = c ( x 1 , x 2 , ⋯   , x n , u ) \sum_{i=1}^{n}b_i(x_1,x_2,\cdots,x_n,u)\cfrac{∂u}{∂x_i}=c(x_1,x_2,\cdots,x_n,u) i=1nbi(x1,x2,,xn,u)xiu=c(x1,x2,,xn,u)

    它的求解归结到 n + 1 n+1 n+1 维一阶线性偏微分方程
    ∑ i = 1 n b i ( x 1 , x 2 , ⋯   , x n , u ) ∂ ϕ ∂ x i + c ( x 1 , x 2 , ⋯   , x n , u ) ∂ ϕ ∂ u = 0 \sum_{i=1}^{n}b_i(x_1,x_2,\cdots,x_n,u)\cfrac{∂\phi}{∂x_i}+c(x_1,x_2,\cdots,x_n,u)\cfrac{∂\phi}{∂u}=0 i=1nbi(x1,x2,,xn,u)xiϕ+c(x1,x2,,xn,u)uϕ=0

    以两自变量 x , y x,y x,y 的一阶拟线性偏微分方程为例
    a u x + b u y = c (3.7) au_x+bu_y=c\tag{3.7} aux+buy=c(3.7)
    其中 a , b , c a,b,c a,b,c 都是变量 x , y , u x,y,u x,y,u 的函数, u = u ( x , y ) u=u(x,y) u=u(x,y) ,且系数 a , b a,b a,b 不同时为零。相应的一阶线性偏微分方程为
    a ∂ ϕ ∂ x + b ∂ ϕ ∂ y + c ∂ ϕ ∂ u = 0 (3.8) a\cfrac{∂\phi}{∂x}+b\cfrac{∂\phi}{∂y}+c\cfrac{∂\phi}{∂u}=0\tag{3.8} axϕ+byϕ+cuϕ=0(3.8)
    ϕ ( x , y , u ) \phi(x,y,u) ϕ(x,y,u) 是方程 (3.8) 的解,则 ϕ ( x , y , u ) = 0 \phi(x,y,u)=0 ϕ(x,y,u)=0 是方程 (3.7) 的隐式解。
    证明: ϕ ( x , y , u ) = 0 \phi(x,y,u)=0 ϕ(x,y,u)=0 给出了隐函数 ϕ ( x , y , u ( x , y ) ) = 0 \phi(x,y,u(x,y))=0 ϕ(x,y,u(x,y))=0 分别对 x , y x,y x,y 求偏导
    ∂ ϕ ∂ x + ∂ ϕ ∂ u ∂ u ∂ x = 0 , ∂ ϕ ∂ y + ∂ ϕ ∂ u ∂ u ∂ y = 0 \cfrac{∂\phi}{∂x}+\cfrac{∂\phi}{∂u}\cfrac{∂u}{∂x}=0,\quad \cfrac{∂\phi}{∂y}+\cfrac{∂\phi}{∂u}\cfrac{∂u}{∂y}=0 xϕ+uϕxu=0,yϕ+uϕyu=0
    从中解出 ∂ ϕ ∂ x = − ∂ ϕ ∂ u ∂ u ∂ x , ∂ ϕ ∂ y = − ∂ ϕ ∂ u ∂ u ∂ y \cfrac{∂\phi}{∂x}=-\cfrac{∂\phi}{∂u}\cfrac{∂u}{∂x},\quad \cfrac{∂\phi}{∂y}=-\cfrac{∂\phi}{∂u}\cfrac{∂u}{∂y} xϕ=uϕxu,yϕ=uϕyu
    带入 (3.8) 可得到方程 (3.7) ,即 ϕ ( x , y , u ) = 0 \phi(x,y,u)=0 ϕ(x,y,u)=0 是方程 (3.7) 的隐式解。反之亦成立。
    一阶线性偏微分方程 (3.8) 的通解可由常微分方程组
    d x a ( x , y , u ) = d y b ( x , y , u ) = d u c ( x , y , u ) (3.9) \cfrac{\mathrm dx}{a(x,y,u)}=\cfrac{\mathrm dy}{b(x,y,u)}=\cfrac{\mathrm du}{c(x,y,u)}\tag{3.9} a(x,y,u)dx=b(x,y,u)dy=c(x,y,u)du(3.9)
    的首次积分确定。上式称为方程 (3.7) 的完全特征方程组,它的积分曲线称为完全特征曲线

    二阶线性偏微分方程的分类和标准式

    以下研究方程的分类,并把各类方程分别化为标准型,这样以后就只需讨论标准型的解法。

    两个自变量方程的分类

    先研究两自变量 x , y x,y x,y 的二阶线性偏微分方程
    a 11 u x x + a 12 u x y + a 22 u y y + b 1 u x + b 2 u y + c u + f = 0 (4.1) a_{11}u_{xx}+a_{12}u_{xy}+a_{22}u_{yy}+b_1u_x+b_2u_y+cu+f=0\tag{4.1} a11uxx+a12uxy+a22uyy+b1ux+b2uy+cu+f=0(4.1)
    其中 a 11 , a 12 , a 22 , b 1 , b 2 , c , f a_{11},a_{12},a_{22},b_1,b_2,c,f a11,a12,a22,b1,b2,c,f 都是自变量 x , y x,y x,y 的函数。我们希望通过适当的自变量变换和未知函数的变换,使方程简化,并在此基础上对方程进行分类。

    作自变量非奇异变换(可逆)
    { ξ = ξ ( x , y ) η = η ( x , y ) \begin{cases} ξ=ξ(x,y) \\ η=η(x,y) \end{cases} {ξ=ξ(x,y)η=η(x,y)
    其中雅可比行列式
    ∣ J ∣ = ∣ ∂ ( ξ , η ) ∂ ( x , y ) ∣ = ∣ ξ x ξ y η x η y ∣ ≠ 0 |J|=|\cfrac{∂(ξ,η)}{∂(x,y)}|=\begin{vmatrix}ξ_x&ξ_y \\ η_x&η_y\end{vmatrix}\neq0 J=(x,y)(ξ,η)=ξxηxξyηy=0
    通过代换可得到新的方程
    A 11 u ξ ξ + 2 A 12 u ξ η + A 22 u η η + B 1 u ξ + B 2 u η + C u + F = 0 (4.2) A_{11}u_{ξξ}+2A_{12}u_{ξη}+A_{22}u_{ηη}+B_1u_ξ+B_2u_η+Cu+F=0\tag{4.2} A11uξξ+2A12uξη+A22uηη+B1uξ+B2uη+Cu+F=0(4.2)
    其中
    { A 11 = a 11 ξ x 2 + 2 a 12 ξ x ξ y + a 22 ξ y 2 A 12 = a 11 ξ x η x + a 12 ( ξ x η y + ξ y η x ) + a 22 ξ y η y A 22 = a 11 η x 2 + 2 a 12 η x η y + a 22 η y 2 B 1 = a 11 ξ x x + 2 a 12 ξ x y + a 22 ξ y y + b 1 ξ x + b 2 ξ y B 2 = a 11 η x x + 2 a 12 η x y + a 22 η y y + b 1 η x + b 2 η y C = c , F = f (4.3) \begin{cases} A_{11}=a_{11}ξ_x^2+2a_{12}ξ_xξ_y+a_{22}ξ_y^2 \\ A_{12}=a_{11}ξ_xη_x+a_{12}(ξ_xη_y+ξ_yη_x)+a_{22}ξ_yη_y \\ A_{22}=a_{11}η_x^2+2a_{12}η_xη_y+a_{22}η_y^2 \\ B_1=a_{11}ξ_{xx}+2a_{12}ξ_{xy}+a_{22}ξ_{yy}+b_1ξ_x+b_2ξ_y \\ B_2=a_{11}η_{xx}+2a_{12}η_{xy}+a_{22}η_{yy}+b_1η_x+b_2η_y \\ C=c,F=f \end{cases} \tag{4.3}