精华内容
下载资源
问答
  • 双曲方程Lax-Wendroff的差分格式程序(Matlab) 一一阶双曲型微分方程的初边值问题 精确解为 二数值解法思想和步骤 2.1:网格剖分 为了用差分方法求解上述问题将求解区域作剖分将空间区间作等分将时间区间作等分并记 ...
  • 偏微分方程数值解 双曲方程-显示与隐式 源代码及算法原理简介 编程语言:Matlab 参考书籍《偏微分方程数值解》
  • 应用matlab处理双曲方程(对流方程)的分数步长法。他的格式见下图。谢谢了![图片](https://img-ask.csdn.net/upload/201605/09/1462799464_3422.jpg)
  • 求解偏微分方程我们常用到差分解法,常见的偏微分方程问题有椭圆型方程第一边值问题、 抛物型方程以及双曲方程等。今天我们学习如何用MATLAB求解偏微分方程。 一维状态空间的偏微分方程MATLAB...
    d24da7f3fe7cd9ea07ccf57aa93300b5.gif

    自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元函数及其偏导数的方程,称之为偏微分方程。

    求解偏微分方程我们常用到差分解法,常见的偏微分方程问题有椭圆型方程第一边值问题、 抛物型方程以及双曲型方程等。 今天我们学习如何用MATLAB求解偏微分方程。

      一维状态空间的偏微分方程的 MATLAB 解法

    工具箱命令介绍

    MATLAB 提供了一个指令 pdepe,用以解以下的 PDE 方程式。

    96adf00f9f28d0c2caa3ac00f3ea47cc.png

    其中时间介于 t0  ≤ t ≤ t f 之间,而位置 x 则介于[a,b]有限区域之间。m 值表示问题的对称性,其可为 0,1 或 2,分别表示平板(slab),圆柱(cylindrical)或球体(spherical)的情形。因而,如果 m > 0 ,则 a 必等于 b ,也就是说其具有圆柱或球体的对称关系。同时,式中 f (x,t,u, ∂u/∂x) 一项为流通量 (flux) ,而 s(x,t,u, ∂u/∂x) 为来源 (source) 项 。c(x,t,u, ∂u/∂x) 为偏微分方程的对角线系数矩阵。若某一对角线元素为 0,则表示该偏微分方程为椭圆型偏微分方程,若为正值(不为 0),则为拋物型偏微分方程。请注意 c 的对角线元素一定不全为 0。偏微分方程初始值可表示为

    6ca9d5ab0639bd21778a9ad09c32c432.png

    而边界条件为

    de4757d0abed12b226d4ac6ebf8e5a3b.png

    其中 x 为两端点位置,即 a 或 b

    用以解含上述初始值及边界值条件的偏微分方程的 MATLAB 命令 pdepe 的用法如下:

    8e506bbdf61409842e68515a0eb8e5cc.png
    • m 为问题之对称参数;

    • xmesh 为空间变量 x 的网格点(mesh)位置向量,即 xmesh = [x0 , x1,L, xN ] ,其中 x0 = a (起点), xN = b (终点)。

    • tspan 为时间变量 t 的向量,即 tspan = [t0 ,t1,L,tM ] ,其中 t0 为起始时间,tM 为终点时间。

    • pdefun 为使用者提供的 pde 函数文件。其函数格式如下:

    565006020183eac15e99d7c672886506.png

    亦即,使用者仅需提供偏微分方程中的系数向量。c , f 和 s 均为行(column)向量,而向量 c 即为矩阵 c 的对角线元素。

    icfun 提供解 u 的起始值,其格式为 u = icfun(x) 值得注意的是 u 为行向量。

    bcfun 使用者提供的边界条件函数,格式如下:

    7eb89336961b122e3947e4570f35d9cf.png

    其中 ul 和 ur 分别表示左边界 (xl = b) 和右边界 (xr = a) u 的近似解。输出变量中, pl和 ql 分别表示左边界 p 和 q 的行向量,而 pr 和 qr 则为右边界 p 和 q 的行向量。

    sol 为解答输出。sol 为多维的输出向量, sol(:,: i) 为 ui 的输出,即 ui  = sol(:,:,i) 。

    元素 ui ( j, k) = sol( j, k,i) 表示在 t = tspan( j) 和 x = xmesh(k) 时 ui 之答案。

    options 为求解器的相关解法参数。详细说明参见 odeset 的使用方法。

    注:

    1. MATLAB PDE 求解器 pdepe 的算法,主要是将原来的椭圆型和拋物线型偏微分方程转化为一组常微分方程。此转换的过程是基于使用者所指定的 mesh 点,以二阶空间离散化(spatial discretization)技术为之(Keel and Berzins,1990),然后以 ode15s 的指令求解。采用 ode15s 的 ode 解法,主要是因为在离散化的过程中,椭圆型偏微分方程被转化为一组代数方程,而拋物线型的偏微分方程则被转化为一组联立的微分方程。因而,原偏微分方程被离散化后,变成一组同时伴有微分方程与代数方程的微分代数方程组,故以 ode15s 便可顺利求解。

    2. x 的取点(mesh)位置对解的精确度影响很大,若 pdepe 求解器给出“…has difficulty finding consistent initial considition”的讯息时,使用者可进一步将 mesh 点取密一点,即增加 mesh 点数。另外,若状态 u 在某些特定点上有较快速的变动时,亦需将此处的点取密集些,以增加精确度。值得注意的是 pdepe 并不会自动做 xmesh 的自动取点,使用者必须观察解的特性,自行作取点的操作。一般而言,所取的点数至少需大于 3 以上。

    3. tspan 的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距(step size)的控制由程序自动完成。

    4. 若要获得特定位置及时间下的解,可配合以 pdeval 命令。使用格式如下:

    9d6a92425f1e539253d412e7612af0d4.png
    • m 代表问题的对称性。m =0 表示平板;m =1 表示圆柱体;m =2 表示球体。其意义同 pdepe 中的自变量 m 。

    • xmesh 为使用者在 pdepe 中所指定的输出点位置向量。xmesh = [ x0 , x1 , L, xN ] 。

    • ui 即 sol( j,:,i) 。也就是说其为 pdepe 输出中第 i 个输出 ui 在各点位置 xmesh 处,时间固定为 t j = tspan( j) 下的解。

    • xout 为所欲内插输出点位置向量。此为使用者重新指定的位置向量。

    • uout 为基于所指定位置 xout ,固定时间 t f 下的相对应输出。

    • duoutdx 为相对应的 dudx 输出值。

    以下将以几个例子,详细说明 pdepe 的用法。

    b228c03480714c8d785430ad7e1570d2.png

     求解一维偏微分方程

    例 1

    试解以下之偏微分方程式

    096f74b3c18821189c7f4eda559fa51d.png

    其中 0 ≤ x ≤ 1 ,且满足以下之条件限制式

    (i)起始值条件

    087bc848f4a9f05e108729d2e8096f5d.png

    注:本问题的解析解为

    e1d63f74533c352b0adce06e771ac7d0.png

    解 

    下面将叙述求解的步骤与过程。当完成以下各步骤后,可进一步将其汇总为一主程序M文件,然后求解。

    步骤 1 

    将欲求解的偏微分方程改写成如式的标准式。

    ea54e785723c2b2a068167d6b0942d03.png80415040e720a5c353d3a832a5e0ab9c.png

    步骤 2

    编写偏微分方程的系数向量函数。

    f18c9f0a30b12dc745861cac792adf8d.pnge2091aee6f073551e4fc25588f1a44c5.png

    步骤 3 

    编写起始值条件。

    f13f1d3d0e0b329e9c0e3c626a469938.png

    步骤 4 

    编写边界条件。在编写之前,先将边界条件改写成标准形式,如式(1), 找出相对应的 p(⋅) 和 q(⋅) 函数,然后写出 MATLAB 的边界条件函数,例如,原边界条件可写成

    fcfc434528f291d372f7cd9762f38c4a.png

    因而,边界条件函数可编写成

    function [pl,ql,pr,qr]=ex20_1bc(xl,ul,xr,ur,t)

    pl=ul;

    ql=0;

    pr=pi*exp(-t);

    qr=1;

    步骤 5

    取点。例如

    x=linspace(0,1,20); %x 取 20 点

    t=linspace(0,2,5); %时间取 5 点输出

    步骤 6

    利用 pdepe 求解。

    m=0; %依步骤 1 之结果

    sol=pdepe(m,@ex20_1pdefun,@ex20_1ic,@ex20_1bc,x,t);

    步骤 7 

    显示结果。

    u=sol(:,:,1);

    surf(x,t,u)

    title('pde 数值解')

    xlabel('位置')

    ylabel('时间' )

    zlabel('u')

    若要显示特定点上的解,可进一步指定 x 或 t 的位置,以便绘图。例如,欲了解时间为 2(终点)时,各位置下的解,可输入以下指令(利用 pdeval 指令):

    figure(2); %绘成图 2

    M=length(t); %取终点时间的下标

    xout=linspace(0,1,100); %输出点位置

    [uout,dudx]=pdeval(m,x,u(M,:),xout);

    plot(xout,uout); %绘图

    title('时间为 2 时,各位置下的解')

    xlabel('x')

    ylabel('u')

    综合以上各步骤,可写成一个程序求解例 1。其参考程序如下:

    向上滑动阅览

    function ex1

    %求解一维热传导偏微分方程的一个综合函数程序

    m=0;

    x=linspace(0,1,20); %xmesh

    t=linspace(0,2,20); %tspan

    %以 pde 求解

    sol=pdepe(m,@ex20_1pdefun,@ex20_1ic,@ex20_1bc,x,t);

    u=sol(:,:,1); %取出答案

    %绘图输出

    figure(1)

    surf(x,t,u)

    title('pde 数值解')

    xlabel('位置 x')

    ylabel('时间 t' )

    zlabel('数值解 u')

    %与解析解做比较

    figure(2)

    surf(x,t,exp(-t)'*sin(pi*x));

    title('解析解')

    xlabel('位置 x')

    ylabel('时间 t' )

    zlabel('数值解 u')

    %t=tf=2 时各位置之解

    figure(3)

    M=length(t); %取终点时间的下表

    xout=linspace(0,1,100); %输出点位置

    [uout,dudx]=pdeval(m,x,u(M,:),xout);

    plot(xout,uout); %绘图

    title('时间为 2 时,各位置下的解')

    xlabel('x')

    ylabel('u')

    %pde 函数

    function [c,f,s]=ex20_1pdefun(x,t,u,dudx)

    c=pi^2;

    f=dudx;

    s=0;

    %初始条件函数

    function u0=ex20_1ic(x)

    u0=sin(pi*x);

    %边界条件函数

    function [pl,ql,pr,qr]=ex20_1bc(xl,ul,xr,ur,t)

    pl=ul;

    ql=0;

    pr=pi*exp(-t);

    qr=1;

    例2

    扩散系统之浓度分布参考如图 1 的装置。管中储放静止液体 B,高度为 L=10 ㎝,放置于充满 A 气体的环境中。假设与B 液体接触面之浓度为 CA0 = 0.01molm3 ,且此浓度不随时间改变而改变,即在操作时间内( h = 10 天)维持定值。气体 A 在液体 B 中之扩散系数为 

    9c42eedf85c7159b68cbeaddded14cff.png

    试决定以下两种情况下,气体 A 溶于液体 B 中之流通量(flux)。

    (a) A 与 B 不发生反应;

    (b) A 与 B 发生以下之反应

    e88d0fef7cf3323c35de56d6f186ef74.png

    其反应速率常数

    62da6b5b57083d0a6b4b772e0bec92ff.png8e62be65aee2668877851abd99603b5c.png

    图1 气体 A 在液体 B 中的扩散

    题意解析:

    (a) 因气体 A 与液体 B 不发生反应,故其扩散现象的质量平衡方程如下:

    575cd4069979f1adb1cbfae80132f33a.png

    依题意,其初始及边界条件为

    4f29f2e21dd473826fe433dbe267d4eb.png

    (b) 在气体 A 与液体 B 会发生一次反应的情况下,其质量平衡方程需改写为

    3f476466a7b2e3bd6988f07681873e21.png

    而起始及边界条件同上。

    在获得浓度分布后,即可以 Fick’s law

    669d889d6d3261e9dee595b11b2f1dec.png

    计算流通量。

    MATLAB 程序设计:

    此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下。

    (a)与标准式96adf00f9f28d0c2caa3ac00f3ea47cc.png比较,可得 C = 1 , f = DAB ∂CA/∂z , s = 0 ,和 m = 0 。

    另外,经与式de4757d0abed12b226d4ac6ebf8e5a3b.png比较后得知,左边界及右边界条件之系数分别为

    左边界( z = 0 ):pl = CA (0,t) − CA0 , ql = 0 。

    右边界( z = L ):pr = 0 , qr =1/DAB

    (b)与标准式96adf00f9f28d0c2caa3ac00f3ea47cc.png比较,可得 m = 0 , C = 1 , f = DAB ∂CA/∂z ,和 s = kCA 。而边界条件之系数同(a)之结果。

    利用以上的处理结果,可编写 MATLAB 参考程序如下:

    向上滑动阅览

    function ex20_3_2

    %扩散系统之浓度分布

    clear

    clc

    global DAB k CA0 

    %给定数据

    CA0=0.01;

    L=0.1; DAB=2e-9; k=2e-7; h=10*24*3600; 

    %取点

    t=linspace(0,h,100);

    z=linspace(0,L,10);

    %case (a) 

    m=0;

    sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t); CA=sol(:,:,1);

    for i=1:length(t) [CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0); NAz(i)=-dCAdz_i*DAB;

    end 

    figure(1) 

    subplot(211) 

    surf(z,t/(24*3600),CA) 

    title('case (a)')

    xlabel('length (m)')

    ylabel('time (day)')

    zlabel('conc. (mol/m^3)')

    subplot(212)

    plot(t/(24*3600),NAz'*24*3600)

    xlabel('time (day)')

    ylabel('flux (mol/m^2.day)')

    %case (b) m=0; sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t); 

    CA=sol(:,:,1);

    for i=1:length(t) [CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0); 

    NAz(i)=-dCAdz_i*DAB;

    end

    %

    figure(2)

    subplot(211)

    surf(z,t/(24*3600),CA) title('case (b)') xlabel('length (m)') ylabel('time (day)') zlabel('conc. (mol/m^3)') subplot(212) plot(t/(24*3600),NAz'*24*3600) xlabel('time (day)') ylabel('flux (mol/m^2.day)') 

    %PDE 函数

    %case (a) 

    function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz) 

    global DAB k CA0

    c=1;

    f=DAB*dCAdz;

    s=0;

    %case (a) function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz) global DAB k CA0

    c=1;

    f=DAB*dCAdz;

    s=k*CA;

    %初始条件函数

    function CA_i=ex20_3_2ic(z)

    CA_i=0; 

    %边界条件函数

    function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t) 

    global DAB k CA0

    pl=CAl-CA0; 

    ql=0; 

    pr=0; 

    qr=1/DAB;

    今天的学习就到这里,下次我们讲讲二维状态空间的偏微分方程的 MATLAB 解法。

      END 

    模友们可能已经发现:现在公众号推送文章的顺序,已经不会按时间排列了。这种变化,可能会让各位模友错过我们每天的推送。

    所以,如果你还想像往常一样,聚焦数模乐园,就需要将“数模乐园”标为星标公众号,同时在阅读完文章后,别忘了给一个“在看”哦。

    星标步骤

    (1)点击页面最上方“数模乐园”,进入公众号主页

    (2)点击右上角的小点点,在弹出页面点击“设为星标”,就可以啦。

    4d4bcb3c7a01b8a8b1927aef5e628159.gif0ed0577648f2ca33e68c2614d4808a31.png

    扫码关注我们

    1f676daabb6440e7c1497cc1a9c3f6f4.png

    2020国际赛QQ交流群

    550eb942c76b8b141b7f8303781016ee.png03a6473a2cdbbc57fd21e9c25909a8f5.gif

    球分享

    03a6473a2cdbbc57fd21e9c25909a8f5.gif

    球点赞

    03a6473a2cdbbc57fd21e9c25909a8f5.gif

    球在看

    展开全文
  • 双曲线差分格式求解,包括一个双曲格式求截断误差及稳定性的例子,以及二维波动方程显格式、交替方向格式ADI1和ADIII格式的求解及matlab 程序,文件中附文档及题目的说明及其解法,易读易懂
  • 求解偏微分方程我们常用到差分解法,常见的偏微分方程问题有椭圆型方程第一边值问题、 抛物型方程以及双曲方程等。今天我们学习如何用MATLAB求解偏微分方程。 一维状态空间的偏微分方程MATLAB...
    bc527590442ea092995a4eb0a329c857.gif

    自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元函数及其偏导数的方程,称之为偏微分方程。

    求解偏微分方程我们常用到差分解法,常见的偏微分方程问题有椭圆型方程第一边值问题、 抛物型方程以及双曲型方程等。 今天我们学习如何用MATLAB求解偏微分方程。

      一维状态空间的偏微分方程的 MATLAB 解法

    工具箱命令介绍

    MATLAB 提供了一个指令 pdepe,用以解以下的 PDE 方程式。

    95c3f710eb74bcb23de3022fcac3ab7d.png

    其中时间介于 t0  ≤ t ≤ t f 之间,而位置 x 则介于[a,b]有限区域之间。m 值表示问题的对称性,其可为 0,1 或 2,分别表示平板(slab),圆柱(cylindrical)或球体(spherical)的情形。因而,如果 m > 0 ,则 a 必等于 b ,也就是说其具有圆柱或球体的对称关系。同时,式中 f (x,t,u, ∂u/∂x) 一项为流通量 (flux) ,而 s(x,t,u, ∂u/∂x) 为来源 (source) 项 。c(x,t,u, ∂u/∂x) 为偏微分方程的对角线系数矩阵。若某一对角线元素为 0,则表示该偏微分方程为椭圆型偏微分方程,若为正值(不为 0),则为拋物型偏微分方程。请注意 c 的对角线元素一定不全为 0。偏微分方程初始值可表示为

    80f6f320dcfdb4bcd97195733691ec9b.png

    而边界条件为

    fcb81fc4d15635d76a712cd1c0246a79.png

    其中 x 为两端点位置,即 a 或 b

    用以解含上述初始值及边界值条件的偏微分方程的 MATLAB 命令 pdepe 的用法如下:

    97ba47ac9963c652c8fce054c76fc67e.png
    • m 为问题之对称参数;

    • xmesh 为空间变量 x 的网格点(mesh)位置向量,即 xmesh = [x0 , x1,L, xN ] ,其中 x0 = a (起点), xN = b (终点)。

    • tspan 为时间变量 t 的向量,即 tspan = [t0 ,t1,L,tM ] ,其中 t0 为起始时间,tM 为终点时间。

    • pdefun 为使用者提供的 pde 函数文件。其函数格式如下:

    70a9130bc27e27a40b24dc2fb7e45c66.png

    亦即,使用者仅需提供偏微分方程中的系数向量。c , f 和 s 均为行(column)向量,而向量 c 即为矩阵 c 的对角线元素。

    icfun 提供解 u 的起始值,其格式为 u = icfun(x) 值得注意的是 u 为行向量。

    bcfun 使用者提供的边界条件函数,格式如下:

    5aebfe635c413984f3ad44390e708476.png

    其中 ul 和 ur 分别表示左边界 (xl = b) 和右边界 (xr = a) u 的近似解。输出变量中, pl和 ql 分别表示左边界 p 和 q 的行向量,而 pr 和 qr 则为右边界 p 和 q 的行向量。

    sol 为解答输出。sol 为多维的输出向量, sol(:,: i) 为 ui 的输出,即 ui  = sol(:,:,i) 。

    元素 ui ( j, k) = sol( j, k,i) 表示在 t = tspan( j) 和 x = xmesh(k) 时 ui 之答案。

    options 为求解器的相关解法参数。详细说明参见 odeset 的使用方法。

    注:

    1. MATLAB PDE 求解器 pdepe 的算法,主要是将原来的椭圆型和拋物线型偏微分方程转化为一组常微分方程。此转换的过程是基于使用者所指定的 mesh 点,以二阶空间离散化(spatial discretization)技术为之(Keel and Berzins,1990),然后以 ode15s 的指令求解。采用 ode15s 的 ode 解法,主要是因为在离散化的过程中,椭圆型偏微分方程被转化为一组代数方程,而拋物线型的偏微分方程则被转化为一组联立的微分方程。因而,原偏微分方程被离散化后,变成一组同时伴有微分方程与代数方程的微分代数方程组,故以 ode15s 便可顺利求解。

    2. x 的取点(mesh)位置对解的精确度影响很大,若 pdepe 求解器给出“…has difficulty finding consistent initial considition”的讯息时,使用者可进一步将 mesh 点取密一点,即增加 mesh 点数。另外,若状态 u 在某些特定点上有较快速的变动时,亦需将此处的点取密集些,以增加精确度。值得注意的是 pdepe 并不会自动做 xmesh 的自动取点,使用者必须观察解的特性,自行作取点的操作。一般而言,所取的点数至少需大于 3 以上。

    3. tspan 的选取主要是基于使用者对那些特定时间的状态有兴趣而选定。而间距(step size)的控制由程序自动完成。

    4. 若要获得特定位置及时间下的解,可配合以 pdeval 命令。使用格式如下:

    b3fd5b3c6fd2ec12424faca251f1433d.png
    • m 代表问题的对称性。m =0 表示平板;m =1 表示圆柱体;m =2 表示球体。其意义同 pdepe 中的自变量 m 。

    • xmesh 为使用者在 pdepe 中所指定的输出点位置向量。xmesh = [ x0 , x1 , L, xN ] 。

    • ui 即 sol( j,:,i) 。也就是说其为 pdepe 输出中第 i 个输出 ui 在各点位置 xmesh 处,时间固定为 t j = tspan( j) 下的解。

    • xout 为所欲内插输出点位置向量。此为使用者重新指定的位置向量。

    • uout 为基于所指定位置 xout ,固定时间 t f 下的相对应输出。

    • duoutdx 为相对应的 dudx 输出值。

    以下将以几个例子,详细说明 pdepe 的用法。

    ca7ffcf3e72480b2a87ccd174c75b3e2.png

     求解一维偏微分方程

    例 1

    试解以下之偏微分方程式

    ca257addae8bcb48d1a83d406211bc64.png

    其中 0 ≤ x ≤ 1 ,且满足以下之条件限制式

    (i)起始值条件

    c45bce8eb0d5cfccfe869ee71100d920.png

    注:本问题的解析解为

    899356c0956f6a36bdf94c3e291b68bf.png

    解 

    下面将叙述求解的步骤与过程。当完成以下各步骤后,可进一步将其汇总为一主程序M文件,然后求解。

    步骤 1 

    将欲求解的偏微分方程改写成如式的标准式。

    c37257f9530a8093b0aca44f37c0dc32.pngefc46151619703e5b933824f9d39a9bd.png

    步骤 2

    编写偏微分方程的系数向量函数。

    cc6e5d845df403d73f1b49910405ec66.png4a1c970753f761cdea473383d9a1271e.png

    步骤 3 

    编写起始值条件。

    feaa2beeced928b11ea878f055353c70.png

    步骤 4 

    编写边界条件。在编写之前,先将边界条件改写成标准形式,如式(1), 找出相对应的 p(⋅) 和 q(⋅) 函数,然后写出 MATLAB 的边界条件函数,例如,原边界条件可写成

    647e19b32845294e2ea717b4b3ed5f4a.png

    因而,边界条件函数可编写成

    function [pl,ql,pr,qr]=ex20_1bc(xl,ul,xr,ur,t)

    pl=ul;

    ql=0;

    pr=pi*exp(-t);

    qr=1;

    步骤 5

    取点。例如

    x=linspace(0,1,20); %x 取 20 点

    t=linspace(0,2,5); %时间取 5 点输出

    步骤 6

    利用 pdepe 求解。

    m=0; %依步骤 1 之结果

    sol=pdepe(m,@ex20_1pdefun,@ex20_1ic,@ex20_1bc,x,t);

    步骤 7 

    显示结果。

    u=sol(:,:,1);

    surf(x,t,u)

    title('pde 数值解')

    xlabel('位置')

    ylabel('时间' )

    zlabel('u')

    若要显示特定点上的解,可进一步指定 x 或 t 的位置,以便绘图。例如,欲了解时间为 2(终点)时,各位置下的解,可输入以下指令(利用 pdeval 指令):

    figure(2); %绘成图 2

    M=length(t); %取终点时间的下标

    xout=linspace(0,1,100); %输出点位置

    [uout,dudx]=pdeval(m,x,u(M,:),xout);

    plot(xout,uout); %绘图

    title('时间为 2 时,各位置下的解')

    xlabel('x')

    ylabel('u')

    综合以上各步骤,可写成一个程序求解例 1。其参考程序如下:

    向上滑动阅览

    function ex1

    %求解一维热传导偏微分方程的一个综合函数程序

    m=0;

    x=linspace(0,1,20); %xmesh

    t=linspace(0,2,20); %tspan

    %以 pde 求解

    sol=pdepe(m,@ex20_1pdefun,@ex20_1ic,@ex20_1bc,x,t);

    u=sol(:,:,1); %取出答案

    %绘图输出

    figure(1)

    surf(x,t,u)

    title('pde 数值解')

    xlabel('位置 x')

    ylabel('时间 t' )

    zlabel('数值解 u')

    %与解析解做比较

    figure(2)

    surf(x,t,exp(-t)'*sin(pi*x));

    title('解析解')

    xlabel('位置 x')

    ylabel('时间 t' )

    zlabel('数值解 u')

    %t=tf=2 时各位置之解

    figure(3)

    M=length(t); %取终点时间的下表

    xout=linspace(0,1,100); %输出点位置

    [uout,dudx]=pdeval(m,x,u(M,:),xout);

    plot(xout,uout); %绘图

    title('时间为 2 时,各位置下的解')

    xlabel('x')

    ylabel('u')

    %pde 函数

    function [c,f,s]=ex20_1pdefun(x,t,u,dudx)

    c=pi^2;

    f=dudx;

    s=0;

    %初始条件函数

    function u0=ex20_1ic(x)

    u0=sin(pi*x);

    %边界条件函数

    function [pl,ql,pr,qr]=ex20_1bc(xl,ul,xr,ur,t)

    pl=ul;

    ql=0;

    pr=pi*exp(-t);

    qr=1;

    例2

    扩散系统之浓度分布参考如图 1 的装置。管中储放静止液体 B,高度为 L=10 ㎝,放置于充满 A 气体的环境中。假设与B 液体接触面之浓度为 CA0 = 0.01molm3 ,且此浓度不随时间改变而改变,即在操作时间内( h = 10 天)维持定值。气体 A 在液体 B 中之扩散系数为 

    79f4e9ba644e4b613fd91ef32601ab16.png

    试决定以下两种情况下,气体 A 溶于液体 B 中之流通量(flux)。

    (a) A 与 B 不发生反应;

    (b) A 与 B 发生以下之反应

    dd4a3e19e191d37b97997d75b4366b5b.png

    其反应速率常数

    e6e43ad0d4b077866e3f96292e79098f.png538044f6ddd94e0b780252b65b240c09.png

    图1 气体 A 在液体 B 中的扩散

    题意解析:

    (a) 因气体 A 与液体 B 不发生反应,故其扩散现象的质量平衡方程如下:

    711c28edcee23c793b428a5065b907da.png

    依题意,其初始及边界条件为

    34d892f6a525b49d40238e0de2b2572d.png

    (b) 在气体 A 与液体 B 会发生一次反应的情况下,其质量平衡方程需改写为

    190bae5201c3550b958e65e281a6a505.png

    而起始及边界条件同上。

    在获得浓度分布后,即可以 Fick’s law

    81f8922927f9df0953e92debf2bed454.png

    计算流通量。

    MATLAB 程序设计:

    此问题依旧可以利用 pdepe 迅速求解。现就各状况的处理过程简述如下。

    (a)与标准式95c3f710eb74bcb23de3022fcac3ab7d.png比较,可得 C = 1 , f = DAB ∂CA/∂z , s = 0 ,和 m = 0 。

    另外,经与式fcb81fc4d15635d76a712cd1c0246a79.png比较后得知,左边界及右边界条件之系数分别为

    左边界( z = 0 ):pl = CA (0,t) − CA0 , ql = 0 。

    右边界( z = L ):pr = 0 , qr =1/DAB

    (b)与标准式95c3f710eb74bcb23de3022fcac3ab7d.png比较,可得 m = 0 , C = 1 , f = DAB ∂CA/∂z ,和 s = kCA 。而边界条件之系数同(a)之结果。

    利用以上的处理结果,可编写 MATLAB 参考程序如下:

    向上滑动阅览

    function ex20_3_2

    %扩散系统之浓度分布

    clear

    clc

    global DAB k CA0 

    %给定数据

    CA0=0.01;

    L=0.1; DAB=2e-9; k=2e-7; h=10*24*3600; 

    %取点

    t=linspace(0,h,100);

    z=linspace(0,L,10);

    %case (a) 

    m=0;

    sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t); CA=sol(:,:,1);

    for i=1:length(t) [CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0); NAz(i)=-dCAdz_i*DAB;

    end 

    figure(1) 

    subplot(211) 

    surf(z,t/(24*3600),CA) 

    title('case (a)')

    xlabel('length (m)')

    ylabel('time (day)')

    zlabel('conc. (mol/m^3)')

    subplot(212)

    plot(t/(24*3600),NAz'*24*3600)

    xlabel('time (day)')

    ylabel('flux (mol/m^2.day)')

    %case (b) m=0; sol=pdepe(m,@ex20_3_2pdefunb,@ex20_3_2ic,@ex20_3_2bc,z,t); 

    CA=sol(:,:,1);

    for i=1:length(t) [CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0); 

    NAz(i)=-dCAdz_i*DAB;

    end

    %

    figure(2)

    subplot(211)

    surf(z,t/(24*3600),CA) title('case (b)') xlabel('length (m)') ylabel('time (day)') zlabel('conc. (mol/m^3)') subplot(212) plot(t/(24*3600),NAz'*24*3600) xlabel('time (day)') ylabel('flux (mol/m^2.day)') 

    %PDE 函数

    %case (a) 

    function [c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz) 

    global DAB k CA0

    c=1;

    f=DAB*dCAdz;

    s=0;

    %case (a) function [c,f,s]=ex20_3_2pdefunb(z,t,CA,dCAdz) global DAB k CA0

    c=1;

    f=DAB*dCAdz;

    s=k*CA;

    %初始条件函数

    function CA_i=ex20_3_2ic(z)

    CA_i=0; 

    %边界条件函数

    function [pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t) 

    global DAB k CA0

    pl=CAl-CA0; 

    ql=0; 

    pr=0; 

    qr=1/DAB;

    今天的学习就到这里,下次我们讲讲二维状态空间的偏微分方程的 MATLAB 解法。

      END 

    模友们可能已经发现:现在公众号推送文章的顺序,已经不会按时间排列了。这种变化,可能会让各位模友错过我们每天的推送。

    所以,如果你还想像往常一样,聚焦数模乐园,就需要将“数模乐园”标为星标公众号,同时在阅读完文章后,别忘了给一个“在看”哦。

    星标步骤

    (1)点击页面最上方“数模乐园”,进入公众号主页

    (2)点击右上角的小点点,在弹出页面点击“设为星标”,就可以啦。

    2ca6852040a8f0ed326c2169289fb031.gifc5bbdeb57ac7360836d5ab435df1c7be.png

    扫码关注我们

    59824b648afe7025559e34707ef4cd5d.png

    2020国际赛QQ交流群

    244e5f868e4b8601e1a71726c1eb35d6.png08a16309b5ec53a381764d98f23acb4b.gif

    球分享

    08a16309b5ec53a381764d98f23acb4b.gif

    球点赞

    08a16309b5ec53a381764d98f23acb4b.gif

    球在看

    展开全文
  • 利用matlab解决这个格式,他的初值条件为在一个区间内,给定一个值。[A,B]区间,a,b为初值。利用matlab做出这个格式。![图片](https://img-ask.csdn.net/upload/201605/09/1462798791_615072.jpg)
  • http://www.mathworks.cn/help/matlab/ref/pdepe.html;jsessionid=143095f6a8df9e5c4b9dfd0f8be0 有很多默认的模式,还不如自己手动解; pdepe Solve initial-boundary value problems for ...

    http://www.mathworks.cn/help/matlab/ref/pdepe.html;jsessionid=143095f6a8df9e5c4b9dfd0f8be0

    有很多默认的模式,还不如自己手动解;


    pdepe

    Solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D

    Syntax

    sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
    sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)
    [sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

    Arguments

    m

    A parameter corresponding to the symmetry of the problem. m can be slab = 0, cylindrical = 1, or spherical = 2.

    pdefun

    A handle to a function that defines the components of the PDE.

    icfun

    A handle to a function that defines the initial conditions.

    bcfun

    A handle to a function that defines the boundary conditions.

    xmesh

    A vector [x0x1, ..., xn] specifying the points at which a numerical solution is requested for every value in tspan. The elements of xmesh must satisfyx0 < x1 < ... < xn. The length of xmesh must be >= 3.

    tspan

    A vector [t0t1, ..., tf] specifying the points at which a solution is requested for every value in xmesh. The elements of tspan must satisfy t0 < t1 < ... < tf. The length of tspan must be >= 3.

    options

    Some options of the underlying ODE solver are available in pdepeRelTolAbsTol,NormControlInitialStepMaxStep, and Events. In most cases, default values for these options provide satisfactory solutions. See odeset for details.

    Description

    sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) solves initial-boundary value problems for systems of parabolic and elliptic PDEs in the one space variable x and time tpdefunicfun, and bcfun are function handles. See the function_handle reference page for more information. The ordinary differential equations (ODEs) resulting from discretization in space are integrated to obtain approximate solutions at times specified in tspan. The pdepe function returns values of the solution on a mesh provided in xmesh.

    Parameterizing Functions explains how to provide additional parameters to the functions pdefunicfun, orbcfun, if necessary.

    pdepe solves PDEs of the form:

    (1-4)

    The PDEs hold for t0 ≤ t ≤ tf and a ≤ x ≤ b. The interval [a,b] must be finite. m can be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry, respectively. If m > 0, then a must be ≥ 0.

    In Equation 1-4f (x,t,u,∂u/∂x) is a flux term and s (x,t,u,∂u/∂x) is a source term. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix c (x,t,u,∂u/∂x). The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if those values of x are mesh points. Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface.

    For t = t0 and all x, the solution components satisfy initial conditions of the form

    (1-5)

    For all t and either x = a or x = b, the solution components satisfy a boundary condition of the form

    (1-6)

    Elements of q are either identically zero or never zero. Note that the boundary conditions are expressed in terms of the flux f rather than ∂u/∂x. Also, of the two coefficients, only p can depend on u.

    In the call sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan):

    • m corresponds to m.

    • xmesh(1) and xmesh(end) correspond to a and b.

    • tspan(1) and tspan(end) correspond to t0 and tf.

    • pdefun computes the terms cf, and s (Equation 1-4). It has the form

      [c,f,s] = pdefun(x,t,u,dudx)
      

      The input arguments are scalars x and t and vectors u and dudx that approximate the solution u and its partial derivative with respect to x, respectively. cf, and s are column vectors. c stores the diagonal elements of the matrix c (Equation 1-4).

    • icfun evaluates the initial conditions. It has the form

      u = icfun(x)
      

      When called with an argument xicfun evaluates and returns the initial values of the solution components at x in the column vector u.

    • bcfun evaluates the terms p and q of the boundary conditions (Equation 1-6). It has the form

      [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
      

      ul is the approximate solution at the left boundary xl = a and ur is the approximate solution at the right boundary xr = bpl and ql are column vectors corresponding to p and q evaluated at xl, similarly prand qr correspond to xr. When m > 0 and a = 0, boundedness of the solution near x = 0 requires that the flux f vanish at a = 0. pdepe imposes this boundary condition automatically and it ignores values returned in pl and ql.

    pdepe returns the solution as a multidimensional array solui = ui = sol(:,:,i) is an approximation to the ith component of the solution vector u. The element ui(j,k) = sol(j,k,i) approximates ui at (t,x) = (tspan(j),xmesh(k)).

    ui = sol(j,:,i) approximates component i of the solution at time tspan(j) and mesh points xmesh(:). Usepdeval to compute the approximation and its partial derivative ∂ui/∂x at points not included in xmesh. Seepdeval for details.

    sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) solves as above with default integration parameters replaced by values in options, an argument created with the odeset function. Only some of the options of the underlying ODE solver are available in pdepeRelTolAbsTolNormControlInitialStep, and MaxStep. The defaults obtained by leaving off the input argument options will generally be satisfactory. See odeset for details.

    [sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) with the 'Events' property inoptions set to a function handle Events, solves as above while also finding where event functionsg(t,u(x,t))are zero. For each function you specify whether the integration is to terminate at a zero and whether the direction of the zero crossing matters. Three column vectors are returned by events:[value,isterminal,direction] = events(m,t,xmesh,umesh)xmesh contains the spatial mesh and umesh is the solution at the mesh points. Use pdeval to evaluate the solution between mesh points. For the I-th event function, value(i) is the value of the function, ISTERMINAL(I) = 1 if the integration is to terminate at a zero of this event function and 0 otherwise. direction(i) = 0 if all zeros are to be computed (the default), +1 if only zeros where the event function is increasing, and -1 if only zeros where the event function is decreasing. Output tsol is a column vector of times specified in tspan, prior to first terminal event. SOL(j,:,:) is the solution at T(j)TE is a vector of times at which events occur. SOLE(j,:,:) is the solution at TE(j) and indices in vector IE specify which event occurred.

    If UI = SOL(j,:,i) approximates component i of the solution at time TSPAN(j) and mesh points XMESHpdevalevaluates the approximation and its partial derivative ∂ui/∂x at the array of points XOUT and returns them inUOUT and DUOUTDX: [UOUT,DUOUTDX] = PDEVAL(M,XMESH,UI,XOUT)

    Note:   The partial derivative ∂ui/∂x is evaluated here rather than the flux. The flux is continuous, but at a material interface the partial derivative may have a jump.

    Examples

    Example 1. This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.

    This equation holds on an interval 0 ≤ x ≤ 1 for times t ≥ 0.

    The PDE satisfies the initial condition

    and boundary conditions

    It is convenient to use local functions to place all the functions required by pdepe in a single function.

    function pdex1
    
    m = 0;
    x = linspace(0,1,20);
    t = linspace(0,2,5);
    
    sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
    % Extract the first solution component as u.
    u = sol(:,:,1);
    
    % A surface plot is often a good way to study a solution.
    surf(x,t,u) 
    title('Numerical solution computed with 20 mesh points.')
    xlabel('Distance x')
    ylabel('Time t')
    
    % A solution profile can also be illuminating.
    figure
    plot(x,u(end,:))
    title('Solution at t = 2')
    xlabel('Distance x')
    ylabel('u(x,2)')
    % --------------------------------------------------------------
    function [c,f,s] = pdex1pde(x,t,u,DuDx)
    c = pi^2;
    f = DuDx;
    s = 0;
    % --------------------------------------------------------------
    function u0 = pdex1ic(x)
    u0 = sin(pi*x);
    % --------------------------------------------------------------
    function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
    pl = ul;
    ql = 0;
    pr = pi * exp(-t);
    qr = 1;
    

    In this example, the PDE, initial condition, and boundary conditions are coded in local functions pdex1pde,pdex1ic, and pdex1bc.

    The surface plot shows the behavior of the solution.

    The following plot shows the solution profile at the final value of t (i.e., t = 2).

    Example 2. This example illustrates the solution of a system of PDEs. The problem has boundary layers at both ends of the interval. The solution changes rapidly for small t.

    The PDEs are

    where F(y) = exp(5.73y) – exp(–11.46y).

    This equation holds on an interval 0 ≤ x ≤ 1 for times t ≥ 0.

    The PDE satisfies the initial conditions

    and boundary conditions

    In the form expected by pdepe, the equations are

    The boundary conditions on the partial derivatives of u have to be written in terms of the flux. In the form expected by pdepe, the left boundary condition is

    and the right boundary condition is

    The solution changes rapidly for small t. The program selects the step size in time to resolve this sharp change, but to see this behavior in the plots, the example must select the output times accordingly. There are boundary layers in the solution at both ends of [0,1], so the example places mesh points near 0 and 1 to resolve these sharp changes. Often some experimentation is needed to select a mesh that reveals the behavior of the solution.

    function pdex4
    m = 0;
    x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
    t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];
    
    sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
    u1 = sol(:,:,1);
    u2 = sol(:,:,2);
    
    figure
    surf(x,t,u1)
    title('u1(x,t)')
    xlabel('Distance x')
    ylabel('Time t')
    
    figure
    surf(x,t,u2)
    title('u2(x,t)')
    xlabel('Distance x')
    ylabel('Time t')
    % --------------------------------------------------------------
    function [c,f,s] = pdex4pde(x,t,u,DuDx)
    c = [1; 1]; 
    f = [0.024; 0.17] .* DuDx; 
    y = u(1) - u(2);
    F = exp(5.73*y)-exp(-11.47*y);
    s = [-F; F]; 
    % --------------------------------------------------------------
    function u0 = pdex4ic(x);
    u0 = [1; 0]; 
    % --------------------------------------------------------------
    function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
    pl = [0; ul(2)]; 
    ql = [1; 0]; 
    pr = [ur(1)-1; 0]; 
    qr = [0; 1]; 
    

    In this example, the PDEs, initial conditions, and boundary conditions are coded in local functions pdex4pde,pdex4ic, and pdex4bc.

    The surface plots show the behavior of the solution components.

    More About

    collapse all

    Tips

    • The arrays xmesh and tspan play different roles in pdepe.

      tspan – The pdepe function performs the time integration with an ODE solver that selects both the time step and formula dynamically. The elements of tspan merely specify where you want answers and the cost depends weakly on the length of tspan.

      xmesh – Second order approximations to the solution are made on the mesh specified in xmesh. Generally, it is best to use closely spaced mesh points where the solution changes rapidly. pdepe doesnot select the mesh in x automatically. You must provide an appropriate fixed mesh in xmesh. The cost depends strongly on the length of xmesh. When m > 0, it is not necessary to use a fine mesh nearx = 0 to account for the coordinate singularity.

    • The time integration is done with ode15spdepe exploits the capabilities of ode15s for solving the differential-algebraic equations that arise when Equation 1-4 contains elliptic equations, and for handling Jacobians with a specified sparsity pattern.

    • After discretization, elliptic equations give rise to algebraic equations. If the elements of the initial conditions vector that correspond to elliptic equations are not "consistent" with the discretization,pdepe tries to adjust them before beginning the time integration. For this reason, the solution returned for the initial time may have a discretization error comparable to that at any other time. If the mesh is sufficiently fine, pdepe can find consistent initial conditions close to the given ones. If pdepedisplays a message that it has difficulty finding consistent initial conditions, try refining the mesh.

      No adjustment is necessary for elements of the initial conditions vector that correspond to parabolic equations.

    References

    [1] Skeel, R. D. and M. Berzins, "A Method for the Spatial Discretization of Parabolic Equations in One Space Variable," SIAM Journal on Scientific and Statistical Computing, Vol. 11, 1990, pp.1–32.

    See Also

    function_handle | ode15s | odeget | odeset | pdeval



    展开全文
  • 4.8.2 偏微分方程在自然科学的很多领域内,都会遇到微分方程初值问题,特别是偏微分方程,它的定解问题是...本小节仅介绍一些最简单、经典的偏微分方程,如椭圆型、双曲型、抛物型等偏微分方程,并给出求解方法。用...

    4.8.2  偏微分方程

    在自然科学的很多领域内,都会遇到微分方程初值问题,特别是偏微分方程,它的定解问题是描述自然界及科学现象的最重要的工具。可以说,几乎自然界和各种现象都可以通过微分方程(特别是偏微分方程)来描述。

    MATLAB提供了一个专门用于求解偏微分方程的工具箱PDE Toolbox。本小节仅介绍一些最简单、经典的偏微分方程,如椭圆型、双曲型、抛物型等偏微分方程,并给出求解方法。用户可以从中了解其解题的基本方法,从而解决类似的问题。

    1.椭圆型问题

    assempde函数是PDE工具箱中的一个基本函数,它使用有限元法组合PDE问题。该函数用来有选择地生成PDE问题的解。可以用assempde函数求解下面的标量椭圆型问题:

    f55dcd6669576c885fbf7d939b977ad8.png

    或系统椭圆型问题:

    005b15c1f7daf840a5699cea7aef59c0.png

    对于标量的情况,用解的列向量代表解矢量u,列矢量中的值对应于p的对应节点处的解。对于具有np个节点的N维系统,u1的前np行描述u的第1个元素,接下来的np行描述u的第2个元素,依次类推。这样,u的元素就作为N块节点行放到u中。assempde函数的调用语法如下。

    (1)u=assempde(b,p,e,t,c,a,f): 通过在线性方程组中剔除Dirichlet边界条件来组合和求解PDE问题。

    (2)[K,F]=assempde(b,p,e,t,c,a,f): 通过刚性位移近似Dirichlet边界条件来组合和求解PDE问题。K和F分别为刚性矩阵和右边项。PDE问题的有限元解为u1=K\F。

    (3)[K,F,B,ud]=assempde(b,p,e,t,c,a,f):通过从线性方程组剔除Dirichlet边界条件来组合PDE问题。u1=K\F,则返回非Dirichlet点上的解。完整的PDE问题可以通过MATLAB中的表达式u=B*u1+ud求解。

    (4)[K,M,F,Q,G,H,R]=assempde(b,p,e,t,c,a,f):给出PDE问题的分离表示。

    (5)u=assempde(K,M,F,Q,G,H,R):将PDE问题的分离表示转换为单一矩阵或矢量的形式,然后通过从线性方程组中剔除Dirichlet边界条件来求解。

    (6)[K1,F1]=assempde(K,M,F,Q,G,H,R):用很大的常数修改Dirichlet边界条件,从而将PDE问题的分离表示转换为单一矩阵或矢量的形式。

    (7)[K1,F1,B,ud]=assempde(K,M,F,Q,G,H,R):从线性方程组中剔除Dirichlet边界条件,从而将PDE问题的分离表示转换为单一矩阵或矢量的形式。

    b描述PDE问题的边界条件。b也可以是边界条件矩阵或边界M文件的文件名。PDE问题几何模型由网络数据p,e,t描述。网格数据的生成可以查询help文档中的initmesh函数。

    系数cadf可以通过多种方式给定。这些系数也可以与时间t相关,在assempde函数中可以看到所有选项的列表。

    4-47  求解L型薄膜的方程a1d57c485362ba2d6d2b864591e8a706.png14f5481026614a4c77fc5bb7506d3df9.pngDirichlet边界条件a5570c34110fa3c559c0fea510034244.png。最后绘图显示结果。

    >> [p,e,t]=initmesh('lshapeg','hmax',0.2);    %生成初始三角形网格,hmax指网格大小

    >> [p,e,t]=refinemesh('lshapeg',p,e,t);       %  将初始的三角形网格细化

    >> u=assempde('lshapeb',p,e,t,1,0,1);         %  求解方程

    >> pdesurf(p,t,u)                                %  绘制结果图形

    lshapeglshapeb分别为表示对象几何模型和边界条件的M文件,为工具箱自带文件。initmesh函数和 refinemesh函数分别对网格模型进行初始化和细化。pdesurf函数绘制解的表面图。L型薄膜泊松方程的解如图4-15所示。

    2.抛物型问题

    MATLAB提供了parabolic函数求解标量抛物型问题:

    508aafc3c3251fce044dd1fdf3758e04.png

    714e720dfd300632e0a04913800671cc.png

    图4-15  L型薄膜泊松方程的解

    或系统PDE问题:

    871388eae4f52d21d202fc5841415e1d.png

    该函数的调用语法如下。

    (1)u1=parabolic(u0,tlist,g,b,p,e,t,c,a,f,d):p,e,t为网格数据,b为边界条件,初值为u0。

    (2)u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d,rtol,atol):atol和rtol为绝对和相对容限。

    (3)u1=parabolic(u0,tlist,K,F,B,ud,M):生成下面PDE问题的解。

    a5ea70bb0f5d2835adc4b4ac5a111d2d.png

    b1334e4beef79fa09ff56320b0f56461.png的初始值为u0

    4-48  求解热传导方程:e08a79d8e7c1fbf455b7debbc0134590.png。其中6159e73b8b9ec5c756b83146f3b368cd.png。在297b77da4dead647f5829aaedff81b40.png的区域内令622d21cc22cbc717f905e73ec94238b4.png,在其他区域内令41ac7d8e3fed74e4d921a6420f1fd96f.png。使用Dirichlet边界条件252c4fe83dbbd5ba9c72ba94c07c2ec4.png。要求计算时间linspace(0,0.1,20)处的解。

    %  生成初始三角形网格数据

    >> [p,e,t]=initmesh('squareg');                %  参数squareg指计算区域是方形的

    >> [p,e,t]=refinemesh('squareg',p,e,t);       %  将初始的三角形网格细化

    >> u0=zeros(size(p,2),1);

    >> ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);    %  查找区域内的元素

    >> u0(ix)=ones(size(ix));                       %  令

    >> tlist=linspace(0,0.1,20);                   %  时间列表

    >> u1=parabolic(u0,tlist,'squareb1',p,e,t,1,0,1,1);   %  求解偏微分方程

    本例首先调用initmesh函数生成偏微分方程的初始网格,然后调用parabolic函数求解偏微分方程,运行结束将显示如下信息:

    96 successful steps

    0 failed attempts

    194 function evaluations

    1 partial derivatives

    20 LU decompositions

    193 solutions of linear systems

    具体的结果u1是一个665*20的矩阵,这里就略去不显示了。

    3.双曲型问题

    MATLAB提供了hyperbolic函数来求解标量双曲型问题:

    6670e3e1b05f3f4a5bbb8c6604e82b94.png

    或系统双曲型问题:

    351cf9fa4b8cea3762bc3bf0076869c6.png

    hyperbolic函数的调用语法如下。

    u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d,rtol,atol)p,e,t为网格数据,b为边界条件,u0为初值,初始导数为ut0atolrtol为绝对和相对容限。

    u1=hyperbolic(u0,ut0,tlist,K,F,B,ud,M)生成下面PDE问题的解。

    edfe3e94ea71d89e17a056fdf1363053.png

    u的初始值为u0ut0

    8bc61dcda49eb17b83c220deea0af044.png

    MATLAB命令窗口输入:

    >> [p,e,t]=initmesh('squareg');                            %   生成初始三角形网格

    >> x=p(1,:)';

    >> y=p(2,:)';

    >> u0=atan(cos(pi/2*x));

    >> ut0=3*sin(pi*x).*exp(cos(pi*y));

    >> tlist=linspace(0,5,31);                                  %  时间列表

    >> uu=hyperbolic(u0,ut0,tlist,'squareb3',p,e,t,1,0,0,1);   %  求解方程

    本例首先调用initmesh函数生成偏微分方程的初始网格,然后调用hyperbolic函数求解偏微分方程,运行结束将显示如下信息:

    462 successful steps

    70 failed attempts

    1066 function evaluations

    1 partial derivatives

    156 LU decompositions

    1065 solutions of linear systems

    具体的结果uu是一个177*31的矩阵,这里就略去不显示了。

    22d0cbf1987404987b9109d9f93bce74.png

    展开全文
  • 本文参考的是来自mooc上北京师范大学彭芳麟老师的计算物理基础基础知识偏微分方程的三种类型椭圆型 初始条件:无抛物型 初始条件:初始温度分布双曲型初始条件:初始位移与初始速度边界条件Dirchlet边界条件区域边界...
  • 本文参考的是来自mooc上北京师范大学彭芳麟老师的计算物理基础基础知识偏微分方程的三种类型椭圆型 初始条件:无抛物型 初始条件:初始温度分布双曲型初始条件:初始位移与初始速度边界条件Dirchlet边界条件区域边界...
  • 介绍了应用最为广泛的椭圆型、双曲型、抛物型偏微分方程的数值解法,而且还详细编程实现了每种方程的多种常见数值解法。 附件使用MATLAB编程来实现这些算法。
  • Hamilton 算子 yj y) Laplace 算子: y2 MATLAB的PDE工具箱能解的方程类型 椭圆型PDE: |(c u) au f in Q为有界平面.c, a, f,和未知量u是标量且是定义在 Q上的 复函数.c可以是 Q上的2X 2的矩阵函数. 抛物型PDE d 十 |...
  • Hamilton算子: Laplace算子: MATLAB的PDE工具箱能解的方程类型: 椭圆型PDE: 为有界平面.c, a, f,和 未知量u 是标量,且是定义在上的 复函数. c 可以是 上的22的矩阵函数. 抛物型PDE 双曲型PDE 特征值问题 d 是 上的复...
  • dsolve 符号计算解微分方程 E e echo M文件被执行指令的显示 edit 启动M文件编辑器 eig 求特征值和特征向量 eigs 求指定的几个特征值 end 控制流FOR等结构体的结尾元素下标 eps 浮点相对精度 error ...
  • 双曲守恒律的ENO格式和WENO格式的MATLAB代码,求解的是Burgers方程的初值问题,分别使用有限体积法4阶ENO格式,有限体积法3阶和5阶WENO方法求解。时间方向用三阶TVD Runge-Kutta方法。研究格式在解光滑情况下和有...
  • 上带有部分Neumann控制和同位观测的双曲方程的适定正则性.证明了由该双曲方程描述的系统不仅在D.Salamon意义下适定的,而且在G.Weiss意义下是正则的,并且证明相关的直接反馈算子为0.本文证明了Rn中三个平行耦合波...

空空如也

空空如也

1 2
收藏数 35
精华内容 14
关键字:

双曲方程matlab

matlab 订阅