关于向量值函数方程变分的一点注记
关于向量值函数方程变分的一点注记
有人说,自己做的都是标量解函数的偏微分方程变分,当遇到解函数是向量值函数时,就不知道有限元空间怎么取,变分也不知道怎么做了。其实是一样的,我这里做个注记。
预备知识
矩阵内积
矩阵内积(Frobenius 含义下的)定义为:
A:B=∑i,jaijbijA: B=\sum_{i, j} a_{i j} b_{i j} A:B=i,j∑aijbij
它表示的是矩阵的对应位置相乘最后再求和。
对于矩阵内积,我们有如下性质,
- DA:B=A:DTBDA:B = A:D^TBDA:B=A:DTB
- AD:B=A:BDTAD:B=A:BD^TAD:B=A:BDT
- A:DB=DTA:BA:DB = D^TA:BA:DB=DTA:B
- A:BD=ADT:BA:BD=AD^T:BA:BD=ADT:B
它从如下表达很容易看出:
A:B=tr(ATB)=tr(ABT)A: B=\operatorname{tr}\left(A^{T} B\right)=\operatorname{tr}\left(A B^{T}\right) A:B=tr(ATB)=tr(ABT)
函数向量的梯度
为了方便和统一,我们以后记标量的梯度是个列向量,而矢量的梯度是个矩阵,每一行都表示矢量的某个分量的梯度,即,形如:
∇f=(fx1fy1fz1fx2fy2fz2)\nabla \mathbf{f}=\left(\begin{array}{ccc} f_{x}^{1} & f_{y}^{1} & f_{z}^{1} \\ f_{x}^{2} & f_{y}^{2} & f_{z}^{2} \end{array}\right) ∇f=(fx1fx2fy1fy2fz1fz2)
有些地方,关于矢量梯度的定义差一个转置。那么,拉普拉斯算子为:
Δf=(Δf1Δf2)\Delta \mathbf{f}=\left(\begin{array}{ccc} \Delta f_{}^{1} \\ \Delta f_{}^{2} \end{array}\right) Δf=(Δf1Δf2)
有些人的这两个的定义跟我这里的定义差一个转置。做问题搞程序的时候,要搞清楚有没有转置。
向量值有限元空间
类似于标量函数的有限元空间,考虑向量值函数有限元空间,以二维向量值为例,可以定义为:
V=span{[φ10],…,[φm0],[0φ1],…,[0φm]}:=span{(ϕi)}\mathbf{V}=\operatorname{span}\left\{\left[\begin{array}{c} \varphi_{1} \\ 0 \end{array}\right], \ldots,\left[\begin{array}{c} \varphi_{m} \\ 0 \end{array}\right],\left[\begin{array}{c} 0 \\ \varphi_{1} \end{array}\right], \ldots,\left[\begin{array}{c} 0 \\ \varphi_{m} \end{array}\right]\right\}:=\operatorname{span}\left\{\left(\phi_{i}\right)\right\} V=span{[φ10],…,[φm0],[0φ1],…,[0φm]}:=span{(ϕi)}
常数矩阵可穿算子
所有用到常矩阵 D\mathbf{D}D 的地方,是可以随意进出梯度和散度算子的,即 ∇Du=D∇u\nabla_{} \mathbf{D} \mathbf{u}=\mathbf{D} \nabla_{} \mathbf{u}∇Du=D∇u ,且 divDu=Ddiv\operatorname{div} \mathbf{D u}=\mathbf{D} \operatorname{div}divDu=Ddiv 。再 一个对于矩阵的散度和梯度算子,是支持分块矩阵的操作的。
举个简单的例子
有了以上的准备,我们举个简单的例子来看向量值解的方程的变分,
我们考虑这样一个例子:
{ut−duuΔu−duvΔv=au(1−u)−buvu+αvt−dvuΔu−dvvΔv=cuvu+α−dv\left\{\begin{array}{l} u_{t}-d_{u u} \Delta_{} u-d_{u v} \Delta_{} v=a u(1-u)-b \frac{u v}{u+\alpha} \\ v_{t}-d_{v u} \Delta_{} u-d_{v v} \Delta_{} v=c \frac{u v}{u+\alpha}-d v \end{array}\right. {ut−duuΔu−duvΔv=au(1−u)−bu+αuvvt−dvuΔu−dvvΔv=cu+αuv−dv
令 u=[u,v]T\mathbf{u}=[u, v]^{T}u=[u,v]T ,那么上面的方程组可以写为:
ut−DΔu=f(u)\mathbf{u}_{t}-D \Delta{} \mathbf{u}=\mathbf{f}(\mathbf{u}) ut−DΔu=f(u)
做变分就得到:
∫ut⋅v+∫∇u:∇(DTv)=∫f(u)⋅v\int \mathbf{u}_{t} \cdot \mathbf{v}+\int \nabla_{} \mathbf{u}: \nabla_{}\left(D^{T} \mathbf{v}\right)=\int \mathbf{f}(\mathbf{u}) \cdot \mathbf{v} ∫ut⋅v+∫∇u:∇(DTv)=∫f(u)⋅v
或者,
∫ut⋅v+∫D∇u:∇v=∫f(u)⋅v\int \mathbf{u}_{t} \cdot \mathbf{v}+\int D\nabla_{} \mathbf{u}: \nabla_{}\mathbf{v}=\int \mathbf{f}(\mathbf{u}) \cdot \mathbf{v} ∫ut⋅v+∫D∇u:∇v=∫f(u)⋅v
PS:
很多人不太理解,为什么 DDD 写在了这个地方,第一个还多了一个转置。这个和我们对向量函数的定义有关。一般来说,我们喜欢把梯度定义为列向量,一列一列地来,但是上述的关于向量值函数的梯度,却是对每个分量求梯度,然后按行排的,这是反直觉的,所以导致了这个地方的 DDD 并不是在 u\mathbf{u}u 前面。
为了和标量函数的梯度是列向量这个事情统一,假设我修改向量值函数梯度和拉普拉斯的定义如下:
∇f=(fx1fy1fz1fx2fy2fz2)T\nabla \mathbf{f}=\left(\begin{array}{ccc} f_{x}^{1} & f_{y}^{1} & f_{z}^{1} \\ f_{x}^{2} & f_{y}^{2} & f_{z}^{2} \end{array}\right)^T ∇f=(fx1fx2fy1fy2fz1fz2)T
Δf=(Δf1Δf2)T\Delta \mathbf{f}=\left(\begin{array}{ccc} \Delta f_{}^{1} \\ \Delta f_{}^{2} \end{array}\right)^T Δf=(Δf1Δf2)T
那么上述问题就应该写为:
ut−(ΔuDT)T=f(u)\mathbf{u}_{t}- (\Delta{} \mathbf{u}D^T)^T=\mathbf{f}(\mathbf{u}) ut−(ΔuDT)T=f(u)
此时我们再做变分,可以得到的是,
∫ut⋅v+∫∇uDT:∇v=∫f(u)⋅v\int \mathbf{u}_{t} \cdot \mathbf{v}+\int \nabla_{} \mathbf{u}D^T: \nabla_{}\mathbf{v}=\int \mathbf{f}(\mathbf{u}) \cdot \mathbf{v} ∫ut⋅v+∫∇uDT:∇v=∫f(u)⋅v
再在 ∇u\nabla_{} \mathbf{u}∇u 和 ∇v\nabla_{} \mathbf{v}∇v 头上加转置变回原来的定义方式,就得到了最开始的上述的表达。
这几行 PS 描述,似乎有点扯淡,这是因为这种情况下,只要保持 ∇u\nabla_{} \mathbf{u}∇u 和 ∇v\nabla_{} \mathbf{v}∇v 定义的方向一致即可,所以不看也罢。
接着考虑向量变分形式的有限元离散,向量有限元空间,定义为:
V=span{[φ10],…,[φm0],[0φ1],…,[0φm]}:=span{(ϕi)}\mathbf{V}=\operatorname{span}\left\{\left[\begin{array}{c} \varphi_{1} \\ 0 \end{array}\right], \ldots,\left[\begin{array}{c} \varphi_{m} \\ 0 \end{array}\right],\left[\begin{array}{c} 0 \\ \varphi_{1} \end{array}\right], \ldots,\left[\begin{array}{c} 0 \\ \varphi_{m} \end{array}\right]\right\}:=\operatorname{span}\left\{\left(\phi_{i}\right)\right\} V=span{[φ10],…,[φm0],[0φ1],…,[0φm]}:=span{(ϕi)}
和标量有限元离散同样的套路,令
u=∑i=1nuiϕi\mathbf{u}=\sum_{i=1}^{n} u_{i} \phi_{i} u=i=1∑nuiϕi
依次取 v=ϕi\mathbf{v}=\phi_{i}v=ϕi ,做向量形式的有限元离散,最后可得:
Mu⃗t+Au⃗=[fu(ϕi)]n×1M \vec{u}_{t}+A \vec{u}=\left[f_{u}\left(\phi_{i}\right)\right]_{n \times 1} Mut+Au=[fu(ϕi)]n×1
其中,
M=[∫ϕi⋅ϕj]n×nM=\left[\int \phi_{i} \cdot \phi_{j}\right]_{n \times n} M=[∫ϕi⋅ϕj]n×n
A=[∫∇ϕi:D∇ϕj]n×nA=\left[\int \nabla_{} \phi_{i}: D\nabla_{} \phi_{j}\right]_{n \times n} A=[∫∇ϕi:D∇ϕj]n×n
fu(ϕi)=∫f(u)⋅ϕif_{u}\left(\phi_{i}\right)=\int \mathbf{f}(\mathbf{u}) \cdot \phi_{i} fu(ϕi)=∫f(u)⋅ϕi
事实上,这个离散系统本质上和对原方程组两个式子分别做变分(用同一个测试空间),最后再拼接组合成一个代数方程组是一样的,即:
(Ms00Ms)u⃗t+(duuAsduvAsdvuAsdvvAs)u⃗=[{f1(u,φi)},{f2(u,φi)}]T\left(\begin{array}{cc} M_{s} & 0 \\ 0 & M_{s} \end{array}\right) \vec{u}_{t}+\left(\begin{array}{cc} d_{u u} A_{s} & d_{u v} A_{s} \\ d_{v u} A_{s} & d_{v v} A_{s} \end{array}\right) \vec{u}=\left[\left\{f_{1}\left(\mathbf{u}, \varphi_{i}\right)\right\},\left\{f_{2}\left(\mathbf{u}, \varphi_{i}\right)\right\}\right]^{T} (Ms00Ms)ut+(duuAsdvuAsduvAsdvvAs)u=[{f1(u,φi)},{f2(u,φi)}]T
这里的 MsM_{s}Ms 和 AsA_{s}As 的下标 sss 表示的是 u,vu, vu,v 所同属的空间 VVV 中的质量矩阵和刚度矩阵。
关于向量值函数方程变分的一点注记相关推荐
- 关于(广义)代数特征值问题的一点注记
关于(广义)代数特征值问题的一点注记 感谢王同学和王同学提供的一些信息和资料. 文章目录 关于(广义)代数特征值问题的一点注记 常用特征值计算方法概述 对称特征值问题 非对称特征值问题 Krylov ...
- CAD插件学习系列教程(十) 高程点注记压盖地物批量处理
大家好,我是王石头,这是CAD插件学习系列教程的第十篇干货文章,这次推荐一款高程点注记压盖地物批量处理插件. 目前使用的数字化成图软件在展绘高程点时,高程注记点存在较多压盖图上已有实体,或者高程注记相 ...
- CAD/CASS高程点自动避让插件、高程点不压线插件(可自动处理高程点注记压线压地物,也可自动处理高程点注记互相压盖)
插件下载:QTools for AutoCAD 功能位置:gcd命令==>高程点避让 功能介绍:可自动处理高程点注记压线压地物,也可自动处理高程点注记互相压盖.(含进度条显示) 功能截图: 视频 ...
- printf格式化输出几点注记
搞了很多年c/c++,有很多细小的东西,曾经不止一次遇到,可是一直都是放在零散的地方,要用的时候怎么也找不到,今天,我痛下决心,改掉不良习惯,把这些经验或是tips记录在这里,便于日后查找. 1.在统 ...
- plot画分段函数_Matlab分段函数图像画法的几点注记
摘要:本文指出<Matlab教程及实训>中关于分段函数画法的不妥,给出Matlab分段函数的最常用的几种画法. 关键词: Matlab 分段函数 Matalb是三大数学软件之一,有强大 ...
- 泛函,变分,欧拉-拉格朗日方程
∫f(Z)p(Z)dZ∫f(Z)p(Z)dZ\int f(Z)p(Z) dZ如何理解?假设Z={z1,z2,...,zn}Z={z1,z2,...,zn}Z=\{z_1,z_2,...,z_n \}. ...
- matlab求异面直线的公垂线,怎么求两异面直线公垂线方程……
共回答了19个问题采纳率:94.7% [思路:先求方向,再求一点.] 记L1和L2的公垂线为L,则L的方程即为所求. 记v为方向向量,记n为法向量. L1:x/1=y/2=z/3 => v(L1 ...
- kcl方程独立性的图论证明
这是困惑了本人很久的问题,最近得闲,去b站观光了部分的离散数学课程并参考了一些资料,对这个问题做个总结,也算是笔记吧.电路原理和离散数学都还没学完,小白一只,各位如有发现任何错误,恳请指出. 支路电流 ...
- matlab求方程在X附近的根,matlab 实验03 求代数方程的近似根(解)
实验三 求代数方程的近似根(解) 求代数方程的根是最常见的数学问题之一(这里称为代数方程,主要是想和后面的微分方程区别开.为简明起见,在本实验的以下叙述中,把代数方程简称为方程),当是一次多项式时,称 ...
最新文章
- 8.4 matlab用户界面设计工具
- html li标签横向排列_HTML简易的常用标签
- 技术分享连载(六十一)
- JEECG开源团队招募新成员 2014年
- VS2005 使用过程中有用的一些快捷键和小技巧
- BUAAOO——UNIT2 SUMMARY
- python爬虫什么意思-网络爬虫是什么(python爬虫有什么用)
- Bejson上线 在线png、jpg图片转svg功能
- 个性签名html,简单个性签名大全2021最新版-简简单单的个性签名--个性简单签名设计-腾牛个性网...
- 每日算法(5)——正整数分解质因数
- 【二】网络空间安全综述
- 大数据学习完整学习路线
- a标签 vue 动态点击_a链接的href怎么用vue动态生成
- The JSP specification requires that an attribute name is preceded by whitespace 解决
- 宝塔linux_开源、强大的Linux服务器集群管理工具,比宝塔好用!
- BMM 和 FMM 中文分词
- 服务器重启后启动php项目
- 5个帮你改掉坏习惯的方法
- VM虚拟机20G磁盘扩展到40G的Linux操作记录
- 基于Android的计步器(Pedometer)的讲解-序