精华内容
下载资源
问答
  • python求解偏微分方程

    2020-04-23 19:33:53
    python求解偏微分方程,画图,容易修改。
  • 基于python偏微分方程自动求解器Fenics说明全书,拥有比tutorial更详细的描述和范例
  • python科学计算解偏微分方程,特别是用高斯迭代法计算拉普拉斯方程 用python科学计算解偏微分方程,特别是用高斯迭代法计算拉普拉斯方程
  • py-pde是一个Python软件包,用于求解偏微分方程(PDE)。 该程序包提供了可在其上定义标量和张量字段的网格的类。 关联的微分算子是使用numba编译的有限差分实现来计算的。 这允许定义,检查和求解典型的PDE,例如...
  • python科学计算解偏微分方程,特别是用高斯迭代法计算拉普拉斯方程
  • Python偏微分方程

    2021-07-14 16:39:57
    泊松问题,包括偏微分方程 − ∇ 2 u = f -\nabla^{2} u=f −∇2u=f 和 ∂ Ω \partial \Omega ∂Ω 上的边界条件 u = u D u=u_{\mathrm{D}} u=uD​,是边界值问题的一个例子,在开始使用 FEniCS 解决它之前必须...

    − ∇ 2 u ( x ) = f ( x ) , x  in  Ω , ( 1 ) -\nabla^{2} u(\boldsymbol{x})=f(\boldsymbol{x}), \quad \boldsymbol{x} \text { in } \Omega,\qquad(1) 2u(x)=f(x),x in Ω,(1)

    u ( x ) = u D ( x ) , x  on  ∂ Ω , ( 2 ) u(\boldsymbol{x})=u_{\mathrm{D}}(\boldsymbol{x}), \quad \boldsymbol{x} \text { on } \partial \Omega,\qquad(2) u(x)=uD(x),x on Ω,(2)

    其中, u = u ( x ) u=u(\boldsymbol{x}) u=u(x) 为未知函数, f = f ( x ) f=f(\boldsymbol{x}) f=f(x) 为规定函数, ∇ 2 \nabla^{2} 2 为拉普拉斯算子(常写为 ∇ \nabla ), Ω \Omega Ω 为空间域, ∂ Ω \partial \Omega Ω Ω \Omega Ω的边界。 泊松问题,包括偏微分方程 − ∇ 2 u = f -\nabla^{2} u=f 2u=f ∂ Ω \partial \Omega Ω 上的边界条件 u = u D u=u_{\mathrm{D}} u=uD,是边界值问题的一个例子,在开始使用 FEniCS 解决它之前必须精确说明。

    在坐标为 x 和 y 的二维空间中,我们可以写出泊松方程为

    − ∂ 2 u ∂ x 2 − ∂ 2 u ∂ y 2 = f ( x , y ) , ( 3 ) -\frac{\partial^{2} u}{\partial x^{2}}-\frac{\partial^{2} u}{\partial y^{2}}=f(x, y),\qquad(3) x22uy22u=f(x,y),(3)

    未知数 u \boldsymbol{u} u 现在是两个变量的函数, u = u ( x , y ) u=u(x, y) u=u(x,y),在二维域 Ω \Omega Ω 上定义。

    泊松方程出现在许多物理环境中,包括热传导、静电、物质扩散、弹性杆的扭曲、无粘性流体流动和水波。 此外,该方程出现在更复杂的偏微分方程系统的数值分裂策略中,尤其是 Navier-Stokes 方程。

    求解边界值问题包括以下步骤:

    有限元变分公式

    FEniCS 基于有限元方法,它是 PDE 数值解的通用且高效的数学机器。 有限元方法的起点是以变分形式表示的偏微分方程。

    将 PDE 转化为变分问题的基本方法是将 PDE 乘以函数 v v v,在域 Ω \Omega Ω 上对所得方程进行积分,并通过具有二阶导数的部分项进行积分。 乘以 PDE 的函数 v v v 称为测试函数。 要逼近的未知函数 u u u 称为试验函数。 术语试验和测试功能也用于 FEniCS 程序。 试验和测试函数属于某些所谓的函数空间,这些函数空间指定了函数的属性

    在本例中,我们首先将泊松方程乘以测试函数 v v v 并在 Ω \Omega Ω 上积分:

    − ∫ Ω ( ∇ 2 u ) v   d x = ∫ Ω f v   d x , ( 4 ) -\int_{\Omega}\left(\nabla^{2} u\right) v \mathrm{~d} x=\int_{\Omega} f v \mathrm{~d} x,\qquad(4) Ω(2u)v dx=Ωfv dx,(4)

    我们在这里让 dx 表示域 Ω 上积分的微分元素。 我们稍后将让 ds 表示在 Ω 边界上积分的微分元素。

    当我们推导出变分公式时,一个常见的规则是我们尽量保持 u u u v v v 的导数的阶数尽可能小。 在这里,我们有 u u u 的二阶空间导数,可以通过应用部分积分技术将其转换为 u u u v v v 的一阶导数。 公式是这样写的

    − ∫ Ω ( ∇ 2 u ) v   d x = ∫ Ω ∇ u ⋅ ∇ v   d x − ∫ ∂ Ω ∂ u ∂ n v   d s , ( 5 ) -\int_{\Omega}\left(\nabla^{2} u\right) v \mathrm{~d} x=\int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x-\int_{\partial \Omega} \frac{\partial u}{\partial n} v \mathrm{~d} s,\qquad(5) Ω(2u)v dx=Ωuv dxΩnuv ds,(5)

    其中 ∂ u ∂ n = ∇ u ⋅ n \frac{\partial u}{\partial n}=\nabla u \cdot n nu=un u u u 在边界外法线方向 n n n 上的导数。

    变分公式的另一个特点是测试函数 v v v 需要在解 u u u 已知的边界部分消失。 在当前问题中,这意味着整个边界 ∂ Ω \partial \Omega Ω 上的 v = 0 v=0 v=0。 因此,等式(5)右边的第二项消失了。 从等式(4)和(5)可以得出

    ∫ Ω ∇ u ⋅ ∇ v   d x = ∫ Ω f v   d x , ( 6 ) \int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x=\int_{\Omega} f v \mathrm{~d} x,\qquad(6) Ωuv dx=Ωfv dx,(6)

    抽象有限元变分公式

    事实证明,为变分问题引入以下规范符号是很方便的:找到 u ∈ V u \in V uV 使得

    a ( u , v ) = L ( v ) ∀ v ∈ V ^ , ( 9 ) a(u, v)=L(v) \quad \forall v \in \hat{V},\qquad(9) a(u,v)=L(v)vV^,(9)

    对于泊松方程,我们有:

    a ( u , v ) = ∫ Ω ∇ u ⋅ ∇ v   d x , ( 10 ) a(u, v)=\int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x,\qquad(10) a(u,v)=Ωuv dx,(10)

    L ( v ) = ∫ Ω f v   d x , ( 11 ) L(v)=\int_{\Omega} f v \mathrm{~d} x,\qquad(11) L(v)=Ωfv dx,(11)

    实现代码(Python)

    from fenics import *
    # Create mesh and define function space
    mesh = UnitSquareMesh(8, 8)
    V = FunctionSpace(mesh, ’P’, 1)
    # Define boundary condition
    u_D = Expression(1 + x[0]*x[0] + 2*x[1]*x[1], degree=2)
    def boundary(x, on_boundary):
    return on_boundary
    bc = DirichletBC(V, u_D, boundary)
    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Constant(-6.0)
    a = dot(grad(u), grad(v))*dx
    L = f*v*dx
    # Compute solution
    u = Function(V)
    solve(a == L, u, bc)
    ...
    

    详情参阅http://viadean.com/py_part_diff_form.html

    展开全文
  • python初探偏微分方程数值解

    千次阅读 2021-01-28 00:32:59
    简要介绍偏微分方程基础概念,之后基于python,对椭圆型偏微分方程数值解进行代码实现

    预备知识

    偏微分方程的三种类型

    • 椭圆形:   − ( ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 ) ϕ ( x , y , z ) = S ( x , y , z ) (1) \ -(\frac{\partial{^2}}{\partial x^2}+\frac{\partial{^2}}{\partial y^2}+\frac{\partial{^2}}{\partial z^2})\phi(x,y,z)=S(x,y,z)\tag{1}  (x22+y22+z22)ϕ(x,y,z)=S(x,y,z)(1)
    • 抛物型:   ( ∂ ∂ t − a 2 ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 ) ϕ ( x , y , z ) = S ( x , y , z ) (2) \ (\frac{\partial}{\partial t}-a^2\frac{\partial{^2}}{\partial x^2}+\frac{\partial{^2}}{\partial y^2}+\frac{\partial{^2}}{\partial z^2})\phi(x,y,z)=S(x,y,z)\tag{2}  (ta2x22+y22+z22)ϕ(x,y,z)=S(x,y,z)(2)
    • 双曲型:   ( ∂ 2 ∂ t 2 − ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 ) ϕ ( x , y , z ) = S ( x , y , z ) (3) \ (\frac{\partial^2}{\partial t^2}-\frac{\partial{^2}}{\partial x^2}+\frac{\partial{^2}}{\partial y^2}+\frac{\partial{^2}}{\partial z^2})\phi(x,y,z)=S(x,y,z)\tag{3}  (t22x22+y22+z22)ϕ(x,y,z)=S(x,y,z)(3)

    定解条件

    偏微分方程的定解条件包括初始条件与边界条件,下面针对三类方程进行简述:

    • 初始条件:
      • 椭圆型:无
      • 抛物型:初始温度分布
      • 双曲型:初始位移与初始速度
    • 边界条件:
      •   D i r i c h l e t \ Dirichlet  Dirichlet边界条件(第一类边界条件):给出区域边界上的函数值。
      •   N e u m a n n \ Neumann  Neumann边界条件(第二类边界条件):给出边界上函数的法向导数。
      • 混合边值条件(第三类边界条件):给出边界上函数及其法向导数的线性组合。

    实战演练

    这里采取对椭圆型方程   u x x + u y y = S ( x , y ) \ u_{xx}+u_{yy}=S(x,y)  uxx+uyy=S(x,y)进行数值求解。

    椭圆型方程离散

      u x x + u y y = S ( x , y ) \ u_{xx}+u_{yy}=S(x,y)  uxx+uyy=S(x,y)进行离散有:
      u ( x + δ x , y ) − 2 u ( x , y ) + u ( x − δ x , y ) ( δ x ) 2 + u ( x , y + δ y ) − 2 u ( x , y ) + u ( x , y − δ y ) ( δ y ) 2 = S ( x , y ) (4) \ \frac{u(x+\delta x,y)-2u(x,y)+u(x-\delta x,y)}{(\delta x)^2}+\frac{u(x,y+\delta y)-2u(x,y)+u(x,y-\delta y)}{(\delta y)^2}=S(x,y)\tag{4}  (δx)2u(x+δx,y)2u(x,y)+u(xδx,y)+(δy)2u(x,y+δy)2u(x,y)+u(x,yδy)=S(x,y)(4)
    其中   δ x = δ y = h \ \delta x=\delta y=h  δx=δy=h
    将其进行整理可得
      u ( i , j ) = 1 4 ( u ( i , j + 1 ) + u ( i , j − 1 ) + u ( i − 1 , j ) + u ( i + 1 , j ) − h 2 S ( i , j ) (5) \ u(i,j)=\frac{1}{4}(u(i,j+1)+u(i,j-1)+u(i-1,j)+u(i+1,j)-h^2S(i,j)\tag{5}  u(i,j)=41(u(i,j+1)+u(i,j1)+u(i1,j)+u(i+1,j)h2S(i,j)(5)
    此时即可采用矩阵法求逆即可,但当网格数过多时,矩阵求逆所需内存量极大,不利于求解,因此引入迭代松弛法进行求解:
      u ( i , j ) = ( 1 − ω ) u ( i , j ) + ω 4 [ u ( i , j + 1 ) + u ( i , j − 1 ) + u ( i − 1 , j ) + u ( i + 1 , j ) − h 2 S ( i , j ) ] (6) \ u(i,j)=(1-\omega)u(i,j)+\\ \frac{\omega}{4}[u(i,j+1)+u(i,j-1)+u(i-1,j)+u(i+1,j)-h^2S(i,j)]\tag{6}  u(i,j)=(1ω)u(i,j)+4ω[u(i,j+1)+u(i,j1)+u(i1,j)+u(i+1,j)h2S(i,j)](6)
    其中   ω \ \omega  ω为松弛因子。
    启动计算时将边界值的平均值作为初始值。初始条件为
      S ( 100 , y ) = 10 , S ( x , 100 ) = 0 , S ( x , 0 ) = 0 , S ( 0 , y ) = 0 \ S(100,y)=10,S(x,100)=0,S(x,0)=0,S(0,y)=0  S(100,y)=10,S(x,100)=0,S(x,0)=0,S(0,y)=0

    python代码

    将上述差分迭代格式以及初始条件代入即可求解,下为python源码:

    '''椭圆形'''
    '''u_xx+u_yy=0'''
    
    h = 0.01
    w = 0.5
    
    u = np.zeros((101,101))
    u[100,:] = 10
    u[1:-1,1:-1] = 2.5
    
    for iteration in range(100):
        for i in range(1,100):
            for j in range(1,100):
                u[i,j] = w/4 * (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1]) + (1-w) * u[i,j]
                
    xl = np.linspace(0, 1, 101)
    yl = np.linspace(0, 1, 101)
    X, Y = np.meshgrid(xl, yl)
    fig = plt.figure()
    ax = Axes3D(fig)
    surf = ax.plot_surface(X, Y, u, cmap=plt.cm.Blues)
    fig.colorbar(surf, shrink=0.5)
    
    

    其中步长设为   h = 0.01 \ h=0.01  h=0.01,松弛因子设为   ω = 0.5 \ \omega=0.5  ω=0.5,其最终结果如图:
    在这里插入图片描述
    在这里插入图片描述可以看出最终的图像是呈凹陷状的,符合客观实际。同时在右侧加入了颜色块便于区分数值。

    以上,请各位巨巨批评指正!

    展开全文
  • 神经网络求解偏微分方程代码分析

    千次阅读 2021-03-20 14:04:15
    下面把我求解偏微分方程的代码分享出来,主要是分享代码思路。这个代码是在求解常微分方程的基础上进行的修改,现在看来有些语句可以换成更高级的表达。若有更好的表达方式欢迎评论/私信! 运行环境:python3.6 + ...

    在我分享了我的神经网络求解微分方程的代码后,很多志同道合的朋友与我进行了交流。下面把我求解偏微分方程的代码分享出来,主要是分享代码思路。这个代码是在求解常微分方程的基础上进行的修改,现在看来有些语句可以换成更高级的表达。若有更好的表达方式欢迎评论/私信!

    运行环境:python3.6 + tensorflow1.2.1 + CPU
    若要tensorflow2.0及以上的版本运行需要添加一行代码

    数学问题

    程序想要解决的数学问题是
    数学公式
    其中求解区域是
    在这里插入图片描述
    对于偏微分方程的 Dirichlet 边界条件由以下方程给出:
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    在这里插入图片描述
    这个微分方程有显式解析解,解析解为:
    在这里插入图片描述
    求解后可以与标准的解析解对比结果。

    代码展示

    import os
    os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'
    #导入相关支持文件包
    import tensorflow as tf
    import matplotlib.pyplot as plt
    import numpy as np
    import math
    from mpl_toolkits.mplot3d import Axes3D
    #产生数据,数据
    x_train_1 = np.linspace(0, 1, 10)#生成[0,1]区间10个点
    x_train_2 = np.linspace(0, 1, 10)#生成[0,1]区间10个点
    input_num = []
    for i in range(len(x_train_1)):
        for j in range(len(x_train_2)):
            input_num.append([x_train_1[i], x_train_2[j]])
    input_num = np.array(input_num)
    input_x1, input_x2 = np.expand_dims(input_num[:, 0], 1), np.expand_dims(input_num[:, 1], 1)
    
    y_trail = np.exp(-1*input_x1)*(input_x1+input_x2**3)
    X2, X1 = np.meshgrid(x_train_1, x_train_2)#使用X1,X2画3维图用
    
    fig = plt.figure("解析解图像")
    ax = fig.gca(projection='3d')
    ax.plot_surface(X1, X2, y_trail.reshape(10, 10), rstride=1, cstride=1, cmap=plt.get_cmap('rainbow'))
    
    def fwd_gradients(Y, x):
        dummy = tf.ones_like(Y)
        G = tf.gradients(Y, x, grad_ys=dummy, colocate_gradients_with_ops=True)[0]
        Y_x = tf.gradients(G, dummy, colocate_gradients_with_ops=True)[0]
        return Y_x
    def act(x):
        return x*tf.nn.sigmoid(x)
    A = tf.placeholder("float", [None, 1])#一次传入100个点[100,1]
    B = tf.placeholder("float", [None, 1])#一次传入100个点[100,1]
    W1_1 = tf.Variable(tf.random_normal([1, 10])*0.01)
    W1_2 = tf.Variable(tf.random_normal([1, 10])*0.01)
    b = tf.Variable(tf.zeros([10])+0.01)
    y1 = act(tf.matmul(A, W1_1)+tf.matmul(B, W1_2)+b)#sigmoid激活函数y1的形状[100,10]
    W2 = tf.Variable(tf.random_normal([10, 1])*0.01)
    y = tf.matmul(y1, W2)#网络的输出[100,1]
    lq = tf.exp(-1*A)*(A-2+B**3+6*B)
    dif_A = fwd_gradients(y, A)
    dif_B = fwd_gradients(y, B)
    dif_AA = fwd_gradients(dif_A, A)
    dif_BB = fwd_gradients(dif_B, B)
    sum = 0
    for i in range(10):
        sum += (y[i, 0]-x_train_2[i]**3)**2
    for i in range(10):
        sum += (y[i+10, 0]-(1+(x_train_2[i])**3)*math.exp(-1))**2
    for i in range(10):
        sum += (y[i*10, 0]-x_train_2[i]*math.exp(-1*x_train_2[i]))**2
    for i in range(10):
        sum += (y[i*10+9, 0]-(x_train_2[i]+1)*math.exp(-1*x_train_2[i]))**2
    t_loss = (dif_AA+dif_BB-lq)**2#常微分方程F的平方
    loss = tf.reduce_mean(t_loss)+sum#未加边界条件!!!!!!!
    train_step = tf.train.AdamOptimizer(0.001).minimize(loss)#Adam优化器训练网络参数
    init = tf.global_variables_initializer()
    with tf.Session() as sess:
        sess.run(init)
        for i in range(50000):#训练50000次
            sess.run(train_step, feed_dict={A: input_x1, B: input_x2})
            if i % 50 == 0:
                total_loss = sess.run(loss, feed_dict={A: input_x1, B: input_x2})
                print("loss={}".format(total_loss))
        output = sess.run(y, feed_dict={A: input_x1, B: input_x2})
        output1 = sess.run(t_loss, feed_dict={A: input_x1, B: input_x2})
        fig = plt.figure("预测曲面")
        bx = fig.gca(projection='3d')
        bx.plot_surface(X1, X2, output.reshape(10, 10), rstride=1, cstride=1, cmap=plt.get_cmap('rainbow'))
        fig = plt.figure("y-y_")
        cset = plt.contourf(X1, X2, (output-y_trail).reshape(10, 10),cmap=plt.get_cmap('rainbow'))
        plt.contour(X1, X2, (output-y_trail).reshape(10, 10))
        plt.colorbar(cset)
        plt.show()
    
    

    代码分析

    程序首先导入相关的支持包文件,然后在[0,1]的区间上产生等距的10个采样点。两个维度分别产生了10个采样点后就能在xoy面上产生100个组合采样点。

    定义的act()方法是激活函数swish。这个函数在求解二阶偏微分方程方面有着较好的性质。除了swish函数,tanh函数也可做为该网络的激活函数。求解偏微分方程涉及到输出对输入的微分,ReLu函数的二阶微分后全都为0,故不能用作激活函数。

    fwd_gradients(Y, x)功能是求出Y对输入x的微分,返回值即为微分后的值。

    结果展示

    在这里插入图片描述
    左图为函数的解析解绘制的图像,右图为神经网络绘制的图像。可以看出通过上述的神经网络能够较好的拟合解析解的结果。

    展开全文
  • 零基础使用 MATLAB 求解偏微分方程(建议收藏) 文章目录零基础使用 MATLAB 求解偏微分方程(建议收藏)偏微分开源工具介绍PDE 工具箱函数汇总介绍0 基础:GUI 界面操作示例问题工具箱求解导出为代码形式代码导出...

    零基础使用 MATLAB 求解偏微分方程(建议收藏)

    偏微分开源工具介绍

    百分之九十以上的重要的工程和数学科学研究,和偏微分方程都脱不开关系。在所有的偏微分方程中,百分之九十九都是没有解析解的。没有解析解怎么办,我们只能通过有限元或者有限差分等方法,求解偏微分方程数值解。如果您有一些代码基础,建议参考我的几篇有限经典博文,简单问题可在此基础上进行修改。

    有限元方法入门:有限元方法简单的一维算例
    有限元方法入门:有限元方法简单的二维算例(三角形剖分)
    有限元方法入门:有限元方法简单的二维算例(矩形剖分)

    对于做工程的朋友,不会偏微分方程数值解,怎么办?没关系,我推荐一些求解各类偏微分方程的容易入门的开源的软件包和工具,它们是:

    • Free FEM++(足够傻瓜又不失自由度,强烈建议做工程的朋友可以学习一下)
    • FEniCS(C++/Python, 开始于芝加哥大学和查尔姆斯理工大学)
    • PETSC (C/Python, 美国阿贡国家实验室)
    • deal.II (C++, 开始于德国海德堡大学)
    • MFEM (C++, 美国劳伦斯利弗莫尔国家实验室)
    • PHG (C, 张林波, 中国科学院)
    • AFEPACK (C++, 李若, 北京大学)
    • FEALPy(Python,魏华祎,湘潭大学)
    • IFEM (MATLAB, 陈龙, UCI)
    • NGSolve(C++/Python,Christoph 等)
    • PHOEBESolver( Fortran and C/C++,宁夏大学,葛永斌)
      [1]: GCGE(C/MATLAB,中国科学院,谢和虎,特征值求解)
    • ……

    当然,商业有限元软件 ansys 等,也非常推荐学工程的朋友去学习,如果要深挖算法的,建议还是用开源的。

    好,有的同学说,这些对你们来说还是太难了。没关系,我可以祭出大招:MATLAB PDE工具箱。为什么它比上面的简单呢?主要是因为,它有可视化的 GUI 工具,你实在不会写代码,你用鼠标点点,也能 “写” 出像模像样的代码。

    PDE 工具箱函数汇总介绍

    PDE 工具箱包含比较多的工具,典型的几个函数如下所示。

    % 偏微分方程工具箱
    %
    % 使用结构化的工作流定义和求解 PDEs
    %   createpde               - 创建 PDE 分析模型
    %   geometryFromEdges       - 从 DECSG or PDEGEOM 创建 2D 几何图形
    %   importGeometry          - 从 STL 文件创建 3D 几何图形
    %   geometryFromMesh        - 从三角网格创建几何图形
    %   multicuboid             - 组合立方体胞单元创建 3D 几何图形
    %   multicylinder           - 组合若干柱状胞单元创建 3D 几何图形
    %   multisphere             - 组合若干球单元创建 3D 几何图形
    %   addVertex               - 在几何区域边界上添加一个顶点
    %   specifyCoefficients     - 指定区域和或者子区域的 PDE 系数
    %   applyBoundaryCondition  - 给几何区域施加边界条件
    %   setInitialConditions    - 设定 PDE 初始条件
    %   generateMesh            - 从几何生成一个网格
    %   solvepde                - 求解 PDE
    %   solvepdeeig             - 求解 PDE 特征值问题
    %   assembleFEMatrices      - 组装中间的有限元矩阵
    %   createPDEResults        - 创建一个用于后处理的结果对象
    %   pdegplot                - 绘制 PDE 几何表示
    %   pdemesh                 - 绘制 PDE 网格
    %   pdeplot                 - 绘制二维 PDE 网格和结果
    %   pdeplot3D               - 绘制 3D PDE 网格和结果
    %   interpolateSolution     - 在指定的空间位置插入解
    %   evaluateGradient        - 在指定的空间位置评估解的梯度
    %   evaluateCGradient       - 评估 PDE 解的通量
    % 
    % 使用热模型解决以传导为主的传热问题
    %   thermalProperties       - 为热模型指定材料的热属性
    %   internalHeatSource      - 指定热模型的内部热源
    %   thermalBC               - 指定热模型的边界条件
    %   thermalIC               - 设置热模型的初始条件或初始猜测
    %   solve                   - 求解热模型中指定的传热问题
    %   interpolateTemperature	- 在任意空间位置的热结果中插入温度
    %   evaluateTemperatureGradient - 评估热解在任意空间位置的温度梯度
    %   evaluateHeatFlux        - 在节点或任意空间位置评估热解的热通量
    %   evaluateHeatRate        - 评估垂直于指定边界的综合热流率
    % 
    % 使用结构模型解决静态、模态和瞬态线性弹性问题
    %   structuralProperties    - 为模型分配结构材料属性
    %   structuralBodyLoad      - 将体载荷应用于结构模型
    %   structuralBoundaryLoad  - 在几何边界上施加结构载荷
    %   structuralBC            - 将边界条件应用于结构模型
    %   structuralIC            - 设置初始位移和速度
    %   structuralDamping       - 为结构模型指定比例阻尼参数
    %   solve                   - 求解 StructuralModel 中指定的结构模型
    %   evaluateStress          - 评估节点位置的应力
    %   evaluateStrain          - E评估节点位置的应变
    %   evaluateVonMisesStress  - 评估节点位置的 von Mises 应力
    %   evaluatePrincipalStrain - 计算节点位置的主要应变
    %   evaluatePrincipalStress - 计算节点位置的主应力
    %   evaluateReaction        - 评估边界上的反作用力
    %   interpolateDisplacement - 在指定的空间位置插入位移
    %   interpolateVelocity     - 在指定的空间位置插入速度
    %   interpolateAcceleration - 在指定的空间位置插入加速度
    %   interpolateStress       - 在指定的空间位置插入应力
    %   interpolateStrain       - 在指定的空间位置插入应变
    %   interpolateVonMisesStress - 在指定空间位置内插 von Mises 应力
    %
    % 使用非结构化工作流程求解 PDE
    %   adaptmesh   - 自适应网格生成和 PDE 解
    %   assema      - 组装面积积分贡献
    %   assemb      - 组装边界条件贡献
    %   assempde    - 组装 PDE 问题
    %   hyperbolic  - 解决双曲线问题
    %   parabolic   - 解决抛物线问题
    %   pdeeig      - 解决特征值 PDE 问题
    %   pdenonlin   - 解决非​​线性 PDE 问题
    %   poisolv     - 矩形网格上泊松方程的快速解
    %
    % 用户界面算法和实用程序
    %   pdecirc     - 画圆
    %   pdeellip    - 绘制椭圆
    %   pdemdlcv    - 转换 MATLAB 4.2c 模型 MATLAB 文件以与 MATLAB 5 一起使用
    %   pdepoly     - 绘制多边形
    %   pderect     - 绘制矩形.
    %   pdeModeler  - PDE Modeler 图形用户界面 (GUI)
    %
    % 几何算法
    %   csgchk      - 检查几何描述矩阵的有效性
    %   csgdel      - 删除最小区域之间的边界
    %   decsg       - 将构造实体几何分解为最小区域
    %   initmesh    - 构建初始三角形网格
    %   jigglemesh  - 抖动三角形网格的内部点
    %   pdearcl     - 参数化表示和弧长之间的插值
    %   poimesh     - 在矩形几何体上制作规则网格
    %   refinemesh  - 细化三角形网格
    %   wbound      - 写入边界条件规范数据文件
    %   wgeom       - W写入几何规格数据文件
    %
    % 绘图函数
    %   pdecont     - 等高线图的速记命令
    %   pdegplot    - 绘制 PDE 几何
    %   pdemesh     - 绘制 PDE 三角形网格
    %   pdeplot     - 通用 PDE 工具箱绘图函数
    %   pdesurf     - 曲面图的速记命令
    %
    % 实用算法
    %   dst         - 离散正弦变换
    %   idst        - 逆离散正弦变换
    %   pdeadgsc    - 使用相对容差标准挑选坏三角形
    %   pdeadworst  - 选择相对于最差值的坏三角形
    %   pdecgrad    - 计算 PDE 解的通量
    %   pdeent      - 与给定三角形集相邻的三角形的索引
    %   pdegrad     - 计算 PDE 解的梯度
    %   pdeintrp    - 将函数值插入到三角形中点
    %   pdejmps     - 适应的误差估计
    %   pdeprtni    - 将函数值内插到网格节点
    %   pdesde      - 与一组子域相邻的边的索引
    %   pdesdp      - 一组子域中的点索引
    %   pdesdt      - 一组子域中的三角形索引
    %   pdesmech    - 计算结构力学张量函数
    %   pdetrg      - 三角形几何数据
    %   pdetriq     - 测量网格三角形的质量
    %   poiasma     - 泊松方程的边界点矩阵贡献
    %   poicalc     - 矩形网格上泊松方程的快速解
    %   poiindex    - 矩形网格的规范排序点的索引
    %   sptarn      - 求解决广义稀疏特征值问题
    %   tri2grid    - 从 PDE 三角形网格插值到矩形网格
    %
    % 用户定义的算法
    %   pdebound    - 边界 MATLAB 文件
    %   pdegeom     - 几何 MATLAB 文件
    % 
    % 对象创建函数。这些函数不是直接调用的。
    %   PDEModel          - 表示 PDE 模型的容器
    %   GeometricModel    - 模型边界的几何表示
    %   AnalyticGeometry  - 来自 PDEGEOM 或 DECSG 几何矩阵的 2D 几何对象
    %   DiscreteGeometry  - 分面边界的几何表示
    %   BoundaryCondition - 定义 PDE 的边界条件
    %   CoefficientAssignmentRecords  - 方程系数的分配
    %   CoefficientAssignment         - 指定区域或子域上的所有 PDE 系数
    %   InitialConditionsRecords   - 记录初始条件的分配
    %   GeometricInitialConditions - 区域或区域边界上的初始条件
    %   NodalInitialConditions  - 在网格节点指定的初始条件
    %   PDEResults           - PDE 解及其派生量
    %   StationaryResults    - PDE 解及其派生量
    %   TimeDependentResults - PDE 解及其派生量
    %   EigenResults         - PDE 解表示
    %   StructuralModel      - 表示结构分析模型的容器
    %   ThermalModel         - 表示热分析模型的容器
    %   ThermalMaterialAssignment - 指定区域或子区域的材料属性
    %   HeatSourceAssignment - 指定域或子域上的热源
    %   ThermalBC            - 定义热模型的边界条件 (BC)
    %   GeometricThermalICs  - 区域或区域边界上的初始温度
    %   NodalThermalICs      - 在网格节点指定的初始温度
    %   ThermalResults       - 热解及其派生量
    %   SteadyStateThermalResults - 稳态热模型解及其派生量
    %   TransientThermalResults   - 瞬态热模型解及其派生量
    %   StructuralMaterialAssignment - 区域或子域上的结构材料属性分配
    %   BodyLoadAssignment   - 结构分析模型的体载荷分配
    %   StructuralBC         - 定义结构模型的边界载荷或边界条件 (BC)
    %   StructuralResults    - 结构解及其派生量
    %   StaticStructuralResults - 静态结构模型解及其派生量
    %   StructuralDampingAssignment - 结构分析模型的阻尼分配
    %   GeometricStructuralICs - 区域上的初始位移和速度
    %   NodalStructuralICs - 在网格节点指定的初始位移和速度
    %   ModalStructuralResults - 结构模态分析结果
    %   TransientStructuralResults - 瞬态结构模型解及其派生量
    %
    
    % 未记录的类和函数
    %   pdeCalcFullU
    %   pdeODEInfo
    %   pdeParabolicInfo
    %   pdeHyperbolicInfo
    
    

    0 基础:GUI 界面操作

    示例问题

    没有什么编程基础,但是又想快速写出有限元程序的同学,建议使用图形界面进行编程,然后导出代码。做个简单的示例操作。比如要求解:
    − Δ u = λ u u ∣ ∂ Ω = 0 Ω 是 一 个 L 型 区 域 , 如 下 图 所 示 -\Delta u = \lambda u\\ u|_{\partial \Omega}=0\\ \Omega 是一个L 型区域,如下图所示 Δu=λuuΩ=0ΩL

    工具箱求解

    • 打开 MATLAB
    • 命令行窗口口输入 pdetool 回车
    • 依次点击菜单栏如下按钮,其中点击 PDE 的时候,改成特征值模式
      在这里插入图片描述

    基础通过 GUI 界面生成的代码此由 pdetool 编写和读取,不应编辑。 有两个推荐的替代方案:

    • 从 pdetool 导出所需的变量并创建一个 MATLAB 脚本,对这些变量执行操作。
    • 使用 MATLAB 脚本完全定义问题。

    导出为代码形式

    得到求解结果后,保存为 main.m 文件,并打开。

    在这里插入图片描述

    function pdemodel
    [pde_fig,ax]=pdeinit;
    pdetool('appl_cb',1);
    set(ax,'DataAspectRatio',[1.5 1 1]);
    set(ax,'PlotBoxAspectRatio',[1 0.74375917767988253 0.74375917767988253]);
    set(ax,'XLim',[-1.5 1.5]);
    set(ax,'YLim',[-1 1]);
    set(ax,'XTickMode','auto');
    set(ax,'YTickMode','auto');
    
    % Geometry description:
    pdepoly([ -1,...
     1,...
     1,...
     0,...
     0,...
     -1,...
    ],...
    [ -1,...
     -1,...
     1,...
     1,...
     0,...
     0,...
    ],...
     'P1');
    set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','P1')
    
    % Boundary conditions:
    pdetool('changemode',0)
    pdesetbd(6,...
    'dir',...
    1,...
    '1',...
    '0')
    pdesetbd(5,...
    'dir',...
    1,...
    '1',...
    '0')
    pdesetbd(4,...
    'dir',...
    1,...
    '1',...
    '0')
    pdesetbd(3,...
    'dir',...
    1,...
    '1',...
    '0')
    pdesetbd(2,...
    'dir',...
    1,...
    '1',...
    '0')
    pdesetbd(1,...
    'dir',...
    1,...
    '1',...
    '0')
    
    % Mesh generation:
    setappdata(pde_fig,'Hgrad',1.3);
    setappdata(pde_fig,'refinemethod','regular');
    setappdata(pde_fig,'jiggle',char('on','mean',''));
    setappdata(pde_fig,'MesherVersion','preR2013a');
    pdetool('initmesh')
    pdetool('refine')
    
    % PDE coefficients:
    pdeseteq(4,...
    '1.0',...
    '0.0',...
    '10.0',...
    '1.0',...
    '0:10',...
    '0.0',...
    '0.0',...
    '[0 100]')
    setappdata(pde_fig,'currparam',...
    ['1.0 ';...
    '0.0 ';...
    '10.0';...
    '1.0 '])
    
    % Solve parameters:
    setappdata(pde_fig,'solveparam',...
    char('0','1548','10','pdeadworst',...
    '0.5','longest','0','1E-4','','fixed','Inf'))
    
    % Plotflags and user data strings:
    setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]);
    setappdata(pde_fig,'colstring','');
    setappdata(pde_fig,'arrowstring','');
    setappdata(pde_fig,'deformstring','');
    setappdata(pde_fig,'heightstring','');
    
    % Solve PDE:
    pdetool('solve') 
    

    代码导出相关数据

    当前目录下,保存如下代码为 matqueque。

    function y = matqueue(p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,p17,p18)
    %MATQUEUE Creates and manipulates a figure-based matrix queue.
    %   FIG = MATQUEUE('create');
    %   Create a queue figure and return its number.
    %
    %   FIG = MATQUEUE('find');
    %   Searches the root window's children to find the queue
    %   figure.  Returns 0 if no queue exists.
    %
    %   MATQUEUE('put', X1, X2, ..., X18);
    %   Insert up to 18 matrices into the queue.  Create the
    %   queue if none exists.
    %   
    %   X = MATQUEUE('get');
    %   Get a matrix out the queue.  Return [] if the queue is
    %   empty.
    %
    %   NUM_ITEMS = MATQUEUE('length');
    %   Return the number of matrices in the queue.  Return -1 if
    %   no buffer exists. 
    %
    %   MATQUEUE('clear')
    %   Empty the queue.
    %
    %   MATQUEUE('close')
    %   Close the queue figure.
    
    
    buffer_name = 'FIFO Buffer';
    
    if (nargin < 1)
      action = 'create';
    else
      action = lower(p0);
    end
    
    if (strcmp(action, 'create'))
      
      %==================================================================
      % Create a new queue.
      %
      % matqueue('create');
      %==================================================================
    
      narginchk(1,1);
      nargoutchk(0,1);
      
      oldFig = findobj(allchild(0), 'flat', 'Visible', 'on');
    
      buffer_fig = matqueue('find');
      if (buffer_fig ~= 0) 
        % Buffer already exists; do nothing
        return;
      end
    
      buffer_fig = figure('Name', buffer_name, ...
          'Visible', 'off',...
          'HandleVisibility', 'callback', ...
          'IntegerHandle', 'off', ...
          'NumberTitle', 'off', ...
          'Tag', buffer_name);
      if (~isempty(oldFig))
        figure(oldFig(1));
      end
      
      queue_holder = uicontrol(buffer_fig, 'Style', 'text', ...
          'Visible', 'off', 'Tag', 'QueueHolder');
    
      y = buffer_fig;
    
      return;
      
    elseif (strcmp(action, 'find'))
      
      %==================================================================
      % Find the queue figure.  If no queue figure exists, return 0.
      %
      % matqueue('find');
      %==================================================================
      
      narginchk(1,1);
      nargoutchk(0,1);
      
      % Search the root's children for a figure with the right name
      buffer_number = findobj(allchild(0), 'flat', 'Tag', buffer_name);
      if (isempty(buffer_number))
        y = 0;
      else
        y = buffer_number(1);
      end
      
      return;
      
    elseif (strcmp(action, 'put'))
      
      %==================================================================
      % Put matrices into the queue.  Queue figure is created if none
      % exists.
      %
      % matqueue('put', X1, X2, ..., X18);
      %==================================================================
      
      narginchk(2,19);
      nargoutchk(0,0);
      
      buffer_fig = matqueue('find');
      if (buffer_fig == 0)
        buffer_fig = matqueue('create');
      end
      
      queue_holder = findobj(get(buffer_fig, 'Children'), 'flat', 'Tag', 'QueueHolder');
      if (isempty(queue_holder))
        error(message('pde:matqueue:corruptMatrixQueue'));
      end
      
      handles = get(queue_holder, 'UserData');
      
      num_inputs = nargin-1;
      new_handles = zeros(1, num_inputs);
      for i = 1:num_inputs
        arg_name = ['p', num2str(i)];
        try_string = ['new_handles(num_inputs+1-i)=uicontrol(buffer_fig,', ...
            ' ''Style'',''text'',''Visible'','...
            ' ''off'',''UserData'',', arg_name, ');'];
        eval(try_string);
      end
      
      set(queue_holder, 'UserData', [new_handles handles]);
      
      return;
      
    elseif (strcmp(action, 'get'))
      
      %==================================================================
      % Return earliest matrix item in the queue.  Errors out if there's
      % no queue.  Returns empty matrix if queue is empty.
      %
      % X = matqueue('get');
      %==================================================================
      
      narginchk(1,1);
      nargoutchk(0,1);
    
      buffer_fig = matqueue('find');
      if (buffer_fig == 0)
        % No buffer; return empty matrix.
        y = [];
        return;
      end
    
      queue_holder = findobj(get(buffer_fig, 'Children'), 'flat', 'Tag', 'QueueHolder');
      if (isempty(queue_holder))
        error(message('pde:matqueue:corruptMatrixQueue'));
      end
      
      handles = get(queue_holder, 'UserData');
      N = length(handles);
      if (N > 0)
        y = get(handles(N), 'UserData');
        delete(handles(N));
        handles(N) = [];
        set(queue_holder, 'UserData', handles);
      else
        % Nothing in the buffer; return empty matrix
        y = [];
      end
      
      return;
      
    elseif (strcmp(action, 'length'))
      
      %==================================================================
      % Returns the length of the queue.  Returns -1 if no queue 
      % figure exists.
      %
      % num_items = matqueue('length');
      %==================================================================
      
      narginchk(1,1);
      nargoutchk(0,1);
      
      buffer_fig = matqueue('find');
      if (buffer_fig == 0)
        % No buffer!  Return 0.
        y = 0;
        return;
      end
      
      queue_holder = findobj(get(buffer_fig, 'Children'), 'flat', 'Tag', 'QueueHolder');
      if (isempty(queue_holder))
        error(message('pde:matqueue:corruptMatrixQueue'));
      end
      
      y = length(get(queue_holder, 'UserData'));
      
      return;
    
    elseif (strcmp(action, 'clear'))
      
      %==================================================================
      % Clear the queue.
      %
      % matqueue('clear');
      %==================================================================
      
      narginchk(1,1);
      nargoutchk(0,1);
      
      buffer_fig = matqueue('find');
      if (buffer_fig == 0)
        % No buffer; nothing to do.
        return;
      end
      
      queue_holder = findobj(get(buffer_fig, 'Children'), 'flat', 'Tag', 'QueueHolder');
      if (~isempty(queue_holder))
        delete(queue_holder);
      end
      
      queue_holder = uicontrol(buffer_fig, 'Style', 'text', ...
          'Visible', 'off', 'Tag', 'QueueHolder');
      
      return;
      
    elseif (strcmp(action, 'close'))
      
      %==================================================================
      % Close the queue figure.
      %
      % matqueue('close');
      %==================================================================
      
      narginchk(1,1);
      nargoutchk(0,1);
      
      buffer_fig = matqueue('find');
      if (buffer_fig ~= 0)
        close(buffer_fig);
      end
      
      return;
      
    else
      
      error(message('pde:matqueue:unknownAction'));
      
    end
    
    

    同样地,将如下代码存为 matqdlg。

    function y = matqdlg(P0,P1,V1,P2,V2,P3,V3,P4,V4,P5,V5,P6,V6,P7,V7,P8,V8,P9,V9)
    %MATQDLG Workspace transfer dialog box.
    %   MATQDLG('ws2buffer', {prop/value pairs})
    %   Put up a dialog box that invites the user to enter a comma-separated
    %   list of expressions.  When user clicks OK,  eval the expressions 
    %   one at a time, putting the results into the buffer.  Allowable
    %   properties include 'PromptString', which may be a string matrix,
    %   'OKCallback', which will be eval'ed when the user finishes
    %   with the dialog box by typing <Return> in the entry field or 
    %   clicking OK, 'CancelCallback', which will be eval'ed when the user 
    %   clicks Cancel, and 'EntryString', the default user entry.  
    %   Figure properties are also allowed in this list, such as 'Name'.
    %
    %   MATQDLG('buffer2ws', {prop/value pairs})
    %   Put up a dialog box that invites the user to enter N comma-separated
    %   variable names, where N is the number of items in the buffer.  Get
    %   items out of the buffer one at a time, storing the results in
    %   indicated workspace variables.  Allowable properties include 
    %   'PromptString', 'OKCallback', 'CancelCallback', and 'EntryString'.
    %   These work the same way as in the 'ws2buffer' action.
    %
    %   Y = MATQDLG('get_entry');
    %   Return the user-entered string in the entry field of the workspace
    %   transfer dialog box.
    %
    %   H = MATQDLG('find')
    %   Return the handle of the dialog box figure.
    %
    %   H = MATQDLG('create')
    %   Create the dialog box figure, returning its handle.
    
    % MATLAB-files required:  matqparse.m, ws2matq.m, matq2ws.m.
    
    buffer_tag = 'Workspace Transfer';
    buffer_name = '';
    
    if (nargin < 1)
      action = 'create';
    else
      action = lower(P0);
    end
    
    if (strcmp(action, 'create'))
      
      %==================================================================
      % Create a new queue.
      %
      % matqdlg('create');
      %==================================================================
    
      narginchk(1,1);
      nargoutchk(0,1);
      
      buffer_fig = matqdlg('find');
      if (buffer_fig ~= 0) 
        % Queue already exists; do nothing
        return;
      end
      
      screenSize = get(0,'ScreenSize');
      width = 415;
      height = 136;
      left = (screenSize(3) - width) / 2;
      bottom = (screenSize(4) - height) / 2;
    
      buffer_fig = figure('Name', buffer_name, 'Visible', 'off',...
          'HandleVisibility', 'callback', ...
          'IntegerHandle', 'off', ...
          'Units', 'pixels', ...
          'Position', [left bottom width height], ...
          'Colormap', [], ...
          'MenuBar', 'none', ...
          'Color', get(0, 'DefaultUIControlBackgroundColor'), ...
          'DefaultUIControlInterruptible','on', ...
          'Tag', buffer_tag, ...
          'NumberTitle', 'off');
      
      axes('Visible', 'off', 'Parent', buffer_fig);
      
      ok = uicontrol(buffer_fig, ...
          'Style', 'push', ...
          'Units', 'normalized', ...
          'Position', [.63 .06 .15 .24], ...
          'Tag', 'OK', ...
          'String', 'OK');
      
      cancel = uicontrol(buffer_fig, ...
          'Style', 'push', ...
          'Units', 'normalized', ...
          'Position', [.80 .06 .15 .24], ...
          'Tag', 'Cancel', ...
          'String', 'Cancel');
      
      prompt = uicontrol(buffer_fig, ...
          'Style', 'text', ...
          'Units', 'normalized', ...
          'Position', [.05 .76 .90 .15], ...
          'String', '', ...
          'Min', 1, ...
          'Max', 3, ...
          'Tag', 'Prompt', ...
          'Horizontal', 'left');
    
      entry = uicontrol(buffer_fig, ...
          'Style', 'edit', ...
          'Units', 'normalized', ...
          'Position', [.05 .46 .90 .31], ...
          'BackgroundColor', 'w', ...
          'ForegroundColor', 'k', ...
          'Tag', 'Entry', ...
          'Horizontal', 'left');
      
      y = buffer_fig;
    
      return;
      
    elseif (strcmp(action, 'find'))
      
      %==================================================================
      % Find the queue figure.  If no queue figure exists, return 0.
      %
      % matqdlg('find');
      %==================================================================
      
      narginchk(1,1);
      nargoutchk(0,1);
      
      % Search the root's children for a figure with the right tag
      buffer_number = findobj(allchild(0), 'flat', 'Type', 'figure', ...
          'Tag', buffer_tag);
      if (isempty(buffer_number))
        y = 0;
      else
        y = buffer_number(1);
      end
      
      return;
      
    elseif (strcmp(action, 'ws2buffer'))
      
      %==================================================================
      % Invoke the dialog box in workspace-to-buffer mode.
      %
      % matqdlg('ws2buffer')
      %==================================================================
      
      narginchk(1,19);
      if (rem(nargin,2) ~= 1)
        error(message('pde:matqdl:invalidNumberOfArgs'));
      end
    
      % Remember the current visible figure.
      figHandles = findobj(allchild(0), 'flat', 'Visible', 'on');
      
      % Set up default properties.
      ok_string = '';
      cancel_string = '';
      entry_string = '';
      prompt_string = 'Enter workspace variable names or expressions:';
      
      buffer_fig = matqdlg('find');
      if (buffer_fig == 0)
        buffer_fig = matqdlg('create');
      end
      if (~isempty(figHandles))
        set(buffer_fig, 'UserData', figHandles(1));
      end
      
      % Process param/value pairs.
      num_properties = (nargin - 1)/2;
      for k = 1:num_properties
        prop_arg_name = ['P' num2str(k)];
        val_arg_name = ['V' num2str(k)];
        prop_arg = lower(eval(prop_arg_name));
        val_arg = eval(val_arg_name);
        if (strcmp(prop_arg, 'promptstring'))
          prompt_string = val_arg;
        elseif (strcmp(prop_arg, 'okcallback'))
          ok_string = val_arg;
        elseif (strcmp(prop_arg, 'cancelcallback'))
          cancel_string = val_arg;
        elseif (strcmp(prop_arg, 'entrystring'))
          entry_string = val_arg;
        else
          set(buffer_fig, prop_arg, val_arg);
        end
      end
      
      promptButton = findobj(buffer_fig, 'Tag', 'Prompt');
      entry = findobj(buffer_fig, 'Tag', 'Entry');
      okButton = findobj(buffer_fig, 'Tag', 'OK');
      cancelButton = findobj(buffer_fig, 'Tag', 'Cancel');
      set(okButton, 'UserData', ok_string, 'Callback', {@okButtonCallback,'ws'});
      set(cancelButton, 'UserData', cancel_string, 'Callback', @cancelButtonCallback);
      set(promptButton, 'String', prompt_string);
      set(entry, 'String', entry_string);
      
      % Adjust height of figure and prompt box.
      prompt_lines = size(prompt_string,1);
      if (prompt_lines > 1)
        set([promptButton entry okButton cancelButton], 'Units', 'pixels');
        prompt_position = get(promptButton, 'Position');
        prompt_height = prompt_position(4);
        increase = (prompt_lines-1) * prompt_height;
        fig_position = get(buffer_fig, 'Position');
        set(promptButton, 'Position', prompt_position + [0 0 0 increase]);
        set(buffer_fig, 'Position', fig_position + [0 0 0 increase]);
        set([promptButton entry okButton cancelButton], 'Units', 'normalized');
      end
    
      drawnow;
      figure(buffer_fig);
      
      return
    
    elseif (strcmp(action, 'buffer2ws'))
      
      %==================================================================
      % Invoke the dialog box in workspace-to-buffer mode.
      %
      % matqdlg('ws2buffer')
      %==================================================================
      
      narginchk(1,19);
      if (rem(nargin,2) ~= 1)
        error(message('pde:matqdl:invalidNumberOfArgs'));
      end
      
      % Remember the current visible figure.
      figHandles = findobj(allchild(0), 'flat', 'Visible', 'on');
      
      num_items = matqueue('length');
      if (num_items == 0)
        % If the queue is empty, there's nothing to do!
        error(message('pde:matqdl:emptyQueue'));
      end
      
      buffer_fig = matqdlg('find');
      if (buffer_fig == 0)
        buffer_fig = matqdlg('create');
      end
      if (~isempty(figHandles))
        set(buffer_fig, 'UserData', figHandles(1));
      end
      
      % Set up default properties.
      ok_string = '';
      entry_string = '';
      cancel_string = '';
      if (num_items == 1)
        prompt_string = 'Enter a variable name:';
      else
        prompt_string = sprintf('Enter %d variable names:', num_items);
      end
      
      % Process param/value pairs.
      num_properties = (nargin - 1)/2;
      for k = 1:num_properties
        prop_arg_name = ['P' num2str(k)];
        val_arg_name = ['V' num2str(k)];
        prop_arg = lower(eval(prop_arg_name));
        val_arg = eval(val_arg_name);
        if (strcmp(prop_arg, 'promptstring'))
          prompt_string = val_arg;
        elseif (strcmp(prop_arg, 'okcallback'))
          ok_string = val_arg;
        elseif (strcmp(prop_arg, 'cancelcallback'))
          cancel_string = val_arg;
        elseif (strcmp(prop_arg, 'entrystring'))
          entry_string = val_arg;
        else
          set(buffer_fig, prop_arg, val_arg);
        end
      end
      
      promptButton = findobj(buffer_fig, 'Tag', 'Prompt');
      set(findobj(buffer_fig, 'Tag', 'OK'), 'UserData', ok_string, ...
          'Callback', {@okButtonCallback,'mat'});
      set(findobj(buffer_fig, 'Tag', 'Cancel'), 'UserData', cancel_string, ...
          'Callback', @cancelButtonCallback);
      set(promptButton, 'String', prompt_string);
      set(findobj(buffer_fig, 'Tag', 'Entry'), 'String', entry_string);
    
      % Adjust height of figure and prompt box.
      prompt_lines = size(prompt_string,1);
      if (prompt_lines > 1)
        set([promptButton entry okButton cancelButton], 'Units', 'pixels');
        prompt_position = get(promptButton, 'Position');
        prompt_height = prompt_position(4);
        increase = (prompt_lines-1) * prompt_height;
        fig_position = get(buffer_fig, 'Position');
        set(promptButton, 'Position', prompt_position + [0 0 0 increase]);
        set(buffer_fig, 'Position', fig_position + [0 0 0 increase]);
        set([promptButton entry okButton cancelButton], 'Units', 'normalized');
      end
    
      drawnow;
      figure(buffer_fig);
      
      return;
      
    
    elseif (strcmp(action, 'get_entry'))
      
      %==================================================================
      % Get the user-entered string in the entry field.
      %
      % string = matqdlg('get_entry');
      %==================================================================
      
      buffer_fig = matqdlg('find');
      if (buffer_fig < 1)
        y = [];
        return;
      end
      
      y = get(findobj(buffer_fig, 'Tag', 'Entry'), 'String');
    
      return;
    
    else
      
      error(message('pde:matqdl:invalidAction'));
      
    end
    
    %--------------------------------------------
    function   okButtonCallback(obj,evd,wsormat)
    
    switch(wsormat)
     case 'ws'
      ws2matq
      
     case 'mat'
      matq2ws
      
     otherwise
      error(message('pde:matqdl:okButtonCallback:unknownOption'));
    
    end
    %-------------------------------------------
    
    function cancelButtonCallback(obj,evd)
    
    matqueue('clear');
    eval(get(findobj(gcbf,'Tag','Cancel'),'UserData'));
    close(matqdlg('find'));
    
    

    将如下代码存为文件 matq2ws 。

    %MATQ2WS Helper script for matqdlg.
    %  MATQ2WS gets the user-entered comma-delimited string,
    %  parses it, and then tries to put the queue contents one
    %  at a time into the resulting variable names.  For
    %  recoverable errors, the prompt is reset and the user can
    %  try again.  Recoverable errors include empty input
    %  string, string containing "#", too few variable names,
    %  and too many variable names.  If the user types something
    %  that cannot be a workspace variable name, that's a
    %  nonrecoverable error.  The queue is cleared and made
    %  invisible, and an error message is printed to the command
    %  window. 
    
    % Variable names (these need to be cleared before returning):
    % var_string_ err_string_ new_prompt_ pound_ N_
    % fatal_error_flag_ i_ expr_ try_string_ catch_string_ error_message_ 
    
    
    var_string_ = get(findobj(gcbf,'Tag','Entry'), 'String');
    [var_string_, err_string_] = matqparse(var_string_);
    if (~isempty(err_string_))
      errordlg(char('Could not parse your expression.', err_string_), ...
          'Workspace Transfer Error', 'on');
      clear var_string_ err_string_ new_prompt_ ...
          pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ ...
          error_message_
      return;
    end
    N_ = size(var_string_, 1);
    if (N_ < matqueue('length'))
      errordlg(char('You did not enter enough variable names.', ...
          'Please try again.'), 'Workspace Transfer Error', 'on');
      clear var_string_ err_string_ new_prompt_ ...
          pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ ...
          error_message_
      return;
    elseif (N_ > matqueue('length'))
      errordlg(char('You entered too many variable names.', ...
          'Please try again.'), 'Workspace Transfer Error', 'on');
      clear var_string_ err_string_ new_prompt_ ...
          pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ ...
          error_message_
      return;
    end
    fatal_error_flag_ = 0;
    for i_ = 1:N_
      expr_ = deblank(var_string_(i_, :));
      try
          assignin('base', expr_, matqueue('get'));
      catch
          fatal_error_flag_ = 1;
      end
    
      if (fatal_error_flag_)
        errordlg(char(sprintf('Error using "%s" as a workspace variable.', ...
      expr_), 'You will need to start over.'), ...
      'Workspace Transfer Error', 'on');
        set(matqdlg('find'), 'Visible', 'off');
        if (~isempty(get(matqdlg('find'),'UserData')))
          if (any(get(0,'Children') == get(matqdlg('find'),'UserData')))
      figure(get(matqdlg('find'),'UserData'));
          end
        end
        matqueue('clear');
        clear var_string_ err_string_ new_prompt_ ...
      pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ ...
      error_message_
        return;
      end
    end
    set(matqdlg('find'), 'Visible', 'off');
    if (~isempty(get(matqdlg('find'),'UserData')))
      if (any(get(0,'Children') == get(matqdlg('find'),'UserData')))
        figure(get(matqdlg('find'),'UserData'));
      end
    end
    
    try
        char(get(get(matqdlg('find'),'CurrentObject'),'UserData'));
    catch
        disp('Error evaluating button callback.')
    end
    close(matqdlg('find'));
    clear var_string_ err_string_ new_prompt_ ...
        pound_ N_ fatal_error_flag_ i_ expr_ try_string_ catch_string_ error_message_ 
    
    

    如下代码为 matqparse 。

    function [m,error_str] = matqparse(str,flag)
    %MATQPARSE Dialog entry parser for MATQDLG.
    %   [M,ERROR_STR] = MATQPARSE(STR,FLAG) is a miniparser
    %   for MATQDLG.
    %   eg: 'abc de  f ghij' becomes [abc ]
    %                                [de  ]
    %                                [f   ]
    %                                [ghij]
    %   Uses either spaces, commas, semi-colons, or brackets
    %   as separators.  Thus 'a 10*[b c] d' will crash. User
    %   must instead say 'a [10*[b c]] d'.
    %
    % See also MATQDLG, MATQUEUE.
    
    % Error checks
    error_str = '';
    if nargin==0
       error_str = getString(message('pde:matqparse:StringReqd'));
       return
    elseif size(str,1)>1 | ~ischar(str)
       error_str = getString(message('pde:matqparse:SingleRowStringReqd'));
       return
    end
    
    if nargin<2
       flag = 1;
    end
    
    l = length(str);
    m = '';
    i = 1;
    j = 1;
    k = 1;
    while k<=l
       % Check for missing [
       if str(k)==']'
          error_str = getString(message('pde:matqparse:UnmatchedRightBracket'));
          return
       elseif str(k)=='['
          % Check for missing ]
          index = find(str(k+1:l)==']');
          if isempty(index)
             error_str = getString(message('pde:matqparse:UnmatchedLeftBracket'));
             return
          else
             % Check for mismatched brackets between k+1 and last element
             index1 = find(str(k+1:l)=='[');
             l_index = length(index);
             l_index1 = length(index1);
             if l_index~=l_index1+1
                error_str = getString(message('pde:matqparse:BracketMismatch'));
                return
             else
                % Everything OK so far
                di = find([index1 index(l_index)+1]>index);
                end_ind = index(di(1));
                m_middle = ['[' matqparse(str(k+1:k+end_ind-1),2) ']'];
                if flag==1
                   % m and m_end may be multiline matrices
                   m_end = matqparse(str(k+end_ind+1:l),1); 
                   m = char(m,m_middle,m_end);
                else
                   % m and m_end will be single line
                   m_end = matqparse(str(k+end_ind+1:l),2);
                   m = [m m_middle m_end];
                end
                k = l+1;
             end
          end
       elseif any(str(k)==' ;,') & (flag==1)
          if j>1
             % Only reset to beginning of next row if 
             % NOT already at beginning of a row
             j=1;
             i = i+1;
          end
          k = k+1; % Increment index into str
       else
          m(i,j) = str(k);
          j = j+1; % Increment column of resultant matrix, m
          k = k+1; % Increment index into str
       end
    end
    
    % Since char of zero is end-of-string flag, change to blanks
    if ~isempty(m)
       EndOfString = find(abs(m)==0);
       m(EndOfString) = char(' '*ones(size(EndOfString)));
    
       % Eliminate any empty rows
       if size(m,2)>1
          m  = m(find(any(m'~=' ')),:);
       end
    
    end
    
    % end matqparse
    
    

    在生成的主文件求解程序末尾添加和修改为如下代码,即可导出数据。

    clc
    clear
    close all
    [pde_fig,ax]=pdeinit;
    pdetool('appl_cb',1);
    set(ax,'DataAspectRatio',[1.5 1 1]);
    set(ax,'PlotBoxAspectRatio',[1 0.74375917767988253 0.74375917767988253]);
    set(ax,'XLim',[-1.5 1.5]);
    set(ax,'YLim',[-1 1]);
    set(ax,'XTickMode','auto');
    set(ax,'YTickMode','auto');
    
    % Geometry description:
    pdepoly([ -1,...
        1,...
        1,...
        0,...
        0,...
        -1,...
        ],...
        [ -1,...
        -1,...
        1,...
        1,...
        0,...
        0,...
        ],...
        'P1');
    set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','P1')
    
    % Boundary conditions:
    pdetool('changemode',0)
    pdesetbd(6,...
        'dir',...
        1,...
        '1',...
        '0')
    pdesetbd(5,...
        'dir',...
        1,...
        '1',...
        '0')
    pdesetbd(4,...
        'dir',...
        1,...
        '1',...
        '0')
    pdesetbd(3,...
        'dir',...
        1,...
        '1',...
        '0')
    pdesetbd(2,...
        'dir',...
        1,...
        '1',...
        '0')
    pdesetbd(1,...
        'dir',...
        1,...
        '1',...
        '0')
    
    % Mesh generation:
    setappdata(pde_fig,'Hgrad',1.3);
    setappdata(pde_fig,'refinemethod','regular');
    setappdata(pde_fig,'jiggle',char('on','mean',''));
    setappdata(pde_fig,'MesherVersion','preR2013a');
    pdetool('initmesh')
    pdetool('refine')
    
    % PDE coefficients:
    pdeseteq(4,...
        '1.0',...
        '0.0',...
        '10.0',...
        '1.0',...
        '0:10',...
        '0.0',...
        '0.0',...
        '[0 100]')
    setappdata(pde_fig,'currparam',...
        ['1.0 ';...
        '0.0 ';...
        '10.0';...
        '1.0 '])
    
    % Solve parameters:
    setappdata(pde_fig,'solveparam',...
        char('0','1548','10','pdeadworst',...
        '0.5','longest','0','1E-4','','fixed','Inf'))
    
    % Plotflags and user data strings:
    setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]);
    setappdata(pde_fig,'colstring','');
    setappdata(pde_fig,'arrowstring','');
    setappdata(pde_fig,'deformstring','');
    setappdata(pde_fig,'heightstring','');
    
    % Solve PDE:
    pdetool('solve')
    for flag=1:6
        % case: export variables to workspace
        pde_fig=findobj(allchild(0),'flat','Tag','PDETool');
        
        if flag==1
            % export geometry data:
            gd=get(findobj(get(pde_fig,'Children'),'flat',...
                'Tag','PDEMeshMenu'),'UserData');
            ns=getappdata(pde_fig,'objnames');
            evalhndl=findobj(get(pde_fig,'Children'),'flat','Tag','PDEEval');
            sf=get(evalhndl,'String');
            matqueue('put',gd,sf,ns)
            pstr='Variable names for geometry data, set formula, labels:';
            estr='gd sf ns';
        elseif flag==2
            % export decomposed list, boundary conditions:
            dl1=getappdata(pde_fig,'dl1');
            h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEBoundMenu');
            bl=get(findobj(get(h,'Children'),'flat',...
                'Tag','PDEBoundMode'),'UserData');
            matqueue('put',dl1,bl)
            pstr='Variable names for decomposed geometry, boundary cond''s:';
            estr='g b';
        elseif flag==3
            % export mesh:
            h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEMeshMenu');
            p=get(findobj(get(h,'Children'),'flat','Tag','PDEInitMesh'),...
                'UserData');
            e=get(findobj(get(h,'Children'),'flat','Tag','PDERefine'),...
                'UserData');
            t=get(findobj(get(h,'Children'),'flat','Tag','PDEMeshParam'),...
                'UserData');
            matqueue('put',p,e,t)
            pstr='Variable names for mesh data (points, edges, triangles):';
            estr='p e t';
        elseif flag==4
            % export PDE coefficients:
            params=get(findobj(get(pde_fig,'Children'),'Tag','PDEPDEMenu'),...
                'UserData');
            ns=getappdata(pde_fig,'ncafd');
            nc=ns(1); na=ns(2); nf=ns(3); nd=ns(4);
            c=params(1:nc,:);
            a=params(nc+1:nc+na,:);
            f=params(nc+na+1:nc+na+nf,:);
            d=params(nc+na+nf+1:nc+na+nf+nd,:);
            matqueue('put',c,a,f,d)
            pstr='Variable names for PDE coefficients:';
            estr='c a f d';
        elseif flag==5
            % export solution:
            u=get(findobj(get(pde_fig,'Children'),'flat','Tag','PDEPlotMenu'),...
                'UserData');
            l=get(findobj(get(pde_fig,'Children'),'flat','Tag','winmenu'),...
                'UserData');
            if isempty(l)
                pstr='Variable name for solution:';
                estr='u';
                matqueue('put',u)
            else
                pstr='Variable names for solution and eigenvalues:';
                estr='u l';
                matqueue('put',u,l)
            end
        elseif flag==6
            % export movie:
            M=getappdata(pde_fig,'movie');
            matqueue('put',M)
            pstr='Variable name for PDE solution movie:';
            estr='M';
        end
        pdeinfo('Change the variable name(s) if desired. OK when done.',0);
        %matqdlg('buffer2ws','Name','Export','PromptString',pstr,...
        %    'OKCallback','pdeinfo;','CancelCallback','pdeinfo;','EntryString',estr);
        
    end
    
    

    出现 You did not enter enough variable names. Please try again. 如何解决?

    不要调用matqdlg('buffer2ws','Name','Export','PromptString',pstr,'OKCallback','pdeinfo;','CancelCallback','pdeinfo;','EntryString',estr);即可。

    0.1 基础:编程调用 PDE 工具箱

    有的朋友说,PDE 工具箱求解 PDE 生成的代码,运行的时候会跳出界面,一看就显得很没有 B 格,有没有办法纯代码操作呢?当然有,界面也只不过是一些代码的集合而已。不想要 PDE 工具箱的界面,又想快速地写出有限元代码,需要一点点代码和有限元基础。

    由工具箱界面生成的代码,一定是和图形界面高度耦合的,我们想把图形界面去掉不显示,并不容易。所以我们考虑使用 MATLAB 脚本完全定义问题。举一个简单的例子来说明它的操作。

    我们还是以拉普拉斯特征值为例:
    − Δ u = λ u -\Delta u=\lambda u Δu=λu

    代码如下:

    clc
    clear
    close all
    model = createpde();%创建PDE模型
    geometryFromEdges(model,@squareg);%从边界生成几何
    pdegplot(model,'EdgeLabels','on')%可视化
    ylim([-1.5,1.5])
    axis equal
    applyBoundaryCondition(model,'dirichlet','Edge',4,'u',0);%左边界 0 狄利克雷边界条件
    applyBoundaryCondition(model,'neumann','Edge',[1,3],'g',0,'q',0);%上下边界 0 纽曼边界条件
    applyBoundaryCondition(model,'neumann','Edge',2,'g',0,'q',-3/4);%由边界混合边界条件,\frac{\partial u}{\partial n}-\frac{3}{4} u=0
    specifyCoefficients(model,'m',0,'d',1,'c',1,'a',0,'f',0);%指定系数,表示特征值问题
    r = [-Inf,10];%找小于10的特征值和特征向量
    generateMesh(model,'Hmax',0.05);%生成网格
    results = solvepdeeig(model,r);%在指定范围求解特征值问题,找六个特征值
    l = results.Eigenvalues;%获得特征值
    u = results.Eigenvectors;%获得特征值对应的特征向量
    pdeplot(model,'XYData',u(:,1));%画第一个特征函数
    pdeplot(model,'XYData',u(:,length(l)));%画最后一个特征函数
    %l(2) - l(1) - pi^2/4
    %l(5) - l(1) - pi^2
    

    上面用的是矩形区域。当然,更复杂的区域我们可以使用 decsg。区域的矩阵表示可以参考这个链接

    decsg 使用的简单示例如下:

    clc
    clear
    close all
    rect1 = [3
        4
        -1
        1
        1
        -1
        0
        0
        -0.5
        -0.5];
    C1 = [1
        1
        -0.25
        0.25];
    C2 = [1
        -1
        -0.25
        0.25];
    C1 = [C1;zeros(length(rect1) - length(C1),1)];
    C2 = [C2;zeros(length(rect1) - length(C2),1)];
    gd = [rect1,C1,C2];
    ns = char('rect1','C1','C2');
    ns = ns';
    sf = '(rect1+C1)-C2';
    [dl,bt] = decsg(gd,sf,ns);
    

    期间的刚度矩阵和质量矩阵也可以通过 assembleFEMatrices 获得。最后给一个区域稍微复杂一点的,通俗易懂的例子吧。顺便和直接把刚度矩阵和质量矩阵导出来调用 eigs 求解做个比较。

    clc
    clear
    close all
    %% 创建模型
    model = createpde;
    radius = 2;
    g = decsg([1 0 0 radius]','C1',('C1')');%通过简单集合图形生成区域
    geometryFromEdges(model,g);
    pdegplot(model,'EdgeLabels','on')%可视化
    axis equal
    title 'Geometry with Edge Labels'
    c = 1;a = 0;f = 0;d = 1;
    specifyCoefficients(model,'m',0,'d',d,'c',c,'a',a,'f',f);
    applyBoundaryCondition(model,'dirichlet','Edge',(1:4),'u',0);
    generateMesh(model,'Hmax',0.2);
    %% 通过导出的刚度矩阵求特征值
    FEMatrices = assembleFEMatrices(model,'nullspace');
    K = FEMatrices.Kc;
    B = FEMatrices.B;
    M = FEMatrices.M;
    sigma = 10; 
    numberEigenvalues = 5;
    [eigenvectorsEigs,eigenvaluesEigs] = eigs(K,M,numberEigenvalues,sigma);
    eigenvaluesEigs = diag(eigenvaluesEigs);
    [maxEigenvaluesEigs,maxIndex] = max(eigenvaluesEigs);
    eigenvectorsEigs = B*eigenvectorsEigs;
    %% 通过工具箱直接求特征值
    r = [min(eigenvaluesEigs)*0.99 max(eigenvaluesEigs)*1.01];
    result = solvepdeeig(model,r);
    eigenvectorsPde = result.Eigenvectors;
    eigenvaluesPde = result.Eigenvalues;
    %% 对比两种方法求出的特征值和特征向量的差距
    eigenValueDiff = sort(eigenvaluesPde) - sort(eigenvaluesEigs);
    fprintf('Maximum difference in eigenvalues from solvepdeeig and eigs: %e\n', ...
      norm(eigenValueDiff,inf));
    %% 画所选范围最大特征值对应的特征函数
    h = figure;
    h.Position = [1 1 2 1].*h.Position;
    subplot(1,2,1)
    axis equal
    pdeplot(model,'XYData', eigenvectorsEigs(:,maxIndex),'Contour','on')
    title(sprintf('eigs eigenvector, eigenvalue: %12.4e', eigenvaluesEigs(maxIndex)))
    xlabel('x')
    ylabel('y')
    subplot(1,2,2)
    axis equal
    pdeplot(model,'XYData',eigenvectorsPde(:,end),'Contour','on')
    title(sprintf('solvepdeeig eigenvector, eigenvalue: %12.4e',eigenvaluesPde(end)))
    xlabel('x')
    ylabel('y')
    
    
    展开全文
  • Python微分方程

    万次阅读 多人点赞 2019-08-19 16:55:45
    Python解微分方程微分方程回顾微分方程:python 解析解(SymPy)微分方程:python数值解(SciPY)微分方程组python数值解 微分方程回顾 微分方程是用来描述某一类函数与其导数之间的关系的方程,其解是一个符合方程...
  • python 求常微分方程 sympy库

    千次阅读 2020-04-26 22:49:34
    利用python的Sympy库求解微分方程的解 y=f(x)y=f(x)y=f(x),并尝试利用matplotlib绘制函数图像 f′(x)+f(x)+f2(x)=0,f(0)=1f'(x)+f(x)+f^2(x)=0,\qquad f(0)=1f′(x)+f(x)+f2(x)=0,f(0)=1 程序,如下 from sympy ...
  • Python小白的数学建模课-11.偏微分方程数值解法

    千次阅读 多人点赞 2021-08-17 14:10:06
    本文采用有限差分法求解偏微分方程,通过案例讲解一维平流方程、一维热传导方程、二维双曲方程、二维抛物方程和二维椭圆方程等常见类型的偏微分方程的数值解法,给出了全部例程和运行结果。。 欢迎关注『Python小白...
  • Galerkin方法求解微分方程组的实现 程序实现
  • 用有限差分法求解一维耦合偏微分方程组 通过极线等值线图显示磁场线的失真 对解的收敛速度和有限差分方法的稳定性进行数值分析 (最终会)演示如何使用 F2PY 极大地提高迭代代码的性能,同时处理 Python 中的所有更...
  • nnde是一组Python模块,这些模块实现了一神经网络,这些神经网络使用技术来求解常微分方程和偏微分方程 。 该软件包在)中进行了描述,该论文的源代码在存储库。 可以在资源库中演示。 接触 埃里克·温特(Eric ...
  • python:使用scipy求解微分方程

    万次阅读 2017-05-06 23:01:36
    遇到一个物理问题,要求解如下微分方程组: d2xdt2=wdydt\frac{d^2x}{dt^2}=w\frac{dy}{dt}d2ydt2=−wdxdt\frac{d^2y}{dt^2}=-w\frac{dx}{dt}经参考相关资料后得知,需要用到scipy包中的odeint函数。 odeint函数...
  • 本文以一个微分方程为例求解微分方程的数值解,并给出模拟数值曲线。普通的微分方程)亦可使用此方法求解并作图。在本文后面部分,给出了python内置函数求解,并比较了二者的区别。
  • matlab优化微分方程组代码图像修复问题的MATLAB / Python代码 这是经典修复方法的详细MATLAB / Python实现。 所提供的所有脚本均在Carola-BibianeSchönlieb所著的“”书中使用,剑桥大学应用数学和计算数学专着,...
  • 数学建模入门-matlab实现偏微分方程数值解

    万次阅读 多人点赞 2019-03-17 22:23:36
    文章目录前言调用示例例题求解命令介绍具体实现步骤1:化标准式步骤 2:编写偏微分方程的系数向量函数步骤3:编写初始条件函数步骤 4:编写边界条件函数步骤 5: 取点主程序 前言 在python3安装fipy失败之后,懒得...
  • COMSOL求解微分方程

    万次阅读 2018-01-10 21:40:16
    COMSOL Multiphysics多物理场仿真软件也提供了求救常微分方程(ODE)和偏微分方程(PDE)的接口,下面详细介绍一下。 (1)建立模型,选择模型向导–>零维–>数学–>全局常微分和微分代数方程(ge),选择研究,选择...
  • 微分方程的解析解(方法归纳)以及基于Python微分方程数值解算例实现 本文归纳常见常微分方程的解析解解法以及基于Python微分方程数值解算例实现。 如果没有出现意外,本文将不包含解法的推导过程。 常微分方程...
  • 通过物理信息神经网络(PINN)求解正向和反向偏微分方程(PDE), 通过PINN求解正向和反向积分微分方程(IDE), 通过分数PINN(fPINN)求解正向和逆分数阶偏微分方程(fPDE), 通过多保真度神经网络(MFNN)从...
  • 背景求解真实问题中建模得到的非线性偏微分方程组, 尽可能少手写代码,如果用matlab, 可能的选项很有限。pdetoolbox有限元方法能够求解一些,但是要求对微分方程的类型以及pdetoolbox特有的书写方式非常熟悉,并能够把...
  • VirtualClaw -- a virtual machine for Clawpack ...Python版PyClaw The 4.6.2 VM contains a minimal Ubuntu distribution together with variouspackages needed for Clawpack, including: gfortran, gcc, g+
  • 椭圆型偏微分方程数值解法

    千次阅读 多人点赞 2020-08-04 19:33:11
    % 一维椭圆方程求解(常微分方程边值问题) % -u'' + q(x)u = f(x), 0<x<1, 取q(x) = x, f(x) = (x-1)exp(x) % u(0) = 1, u(1) = e; 边界条件 % 真解为 u = exp(x) N = 20; h = 1/N; x_all = (0:h:1)'; x = x_...
  • 神经网络求解二阶常微分方程

    千次阅读 热门讨论 2020-11-20 15:10:41
    需要让我使用深度神经网络求解偏微分方程。在相关调研过程中,CSDN上作者Trytobenice分享过相关的程序源码。基于相关程序源码,我将他的一阶常微分方程求解扩充到二阶常微分方程求解。并且按照此方法可以求解高阶常...
  • 数学建模-用Sympy库求解一些微分方程 例题1:解下列微分方程的特解 from sympy.abc import x #因为sympy的计算都是需要符号变量声明的 例如x=symbols('x'); 使用abc库可以将任何拉丁和希腊字母导出为符号型的 from ...
  • 用RK45解偏微分方程,变步长,注释详细。
  • scipy求解微分方程(scipy, torchdiffeq)

    万次阅读 多人点赞 2018-06-13 00:51:52
    使用scipy求解微分方程主要使用scipy.integrate模块,函数是odeint,solve_ivp(初值问题),可以求解一阶、二阶以及高阶方程或方程组。 下面直接上代码,已有详细注释 ''' 使用scipy求解微分方程,包括...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 881
精华内容 352
关键字:

python求解偏微分方程组

python 订阅