Mathematica 之微分特征系统

1 1

导读

版本 11 扩展了其符号和数值微分方程的求解功能，其中包括了在一定区域上寻找特征值和特征函数. 给定一个有可能耦合的偏微分方程（PDE），一个指定的区域，或者再加上边界条件，特征值求解器将求解该 PDE 算符在给定区域上相应的特征值和特征函数. 此功能将使你对物理和工程设计有更深入的了解。

2 1

案例

Mathematica关于微分特征系统的应用部分示例如下:

下面小编用Mathematica求解几个实例的过程向大家展示其在微分特征系统中的应用.

示例1:模拟 CO 分子上的小振动

从实验上看，一个 CO 分子会以有效弹性系数

在其平衡长度附近振动. 振动由量子简谐振子方程描述. 在下面的例子中，

是分子的折合质量，

为固有频率，

为相对平衡位置的位移，

为约化普朗克常数.

In[1]:=qho = -(\[HBar]^2/(2 m)) Laplacian[u[x], {x}] + (m \[Omega]^2)/2 x^2 u[x];

计算前四个特征值和归一化特征函数.

In[2]:=sol = DEigensystem[qho, u[x], {x, -\[Infinity], \[Infinity]},4, Assumptions->\[HBar] > 0 && m > 0 && \[Omega] > 0,Method ->"Normalize"]

假定粒子处于四个态的等量叠加，则波函数将有

的形式.

In[3]:=\[Psi][x_, t_] = Total[MapThread[1/2 Exp[I E t #1/\[HBar]] #2 &, sol]]

用原子质量单位、飞秒和皮米等基本单位来计算三个常数

、

和

，这样计算的值将接近 1 的量级.

In[4]:=m = QuantityMagnitude[(Entity["Element", "Carbon"][EntityProperty["Element", "AtomicMass"]] Entity["Element","Oxygen"][EntityProperty["Element", "AtomicMass"]])/(Entity["Element", "Carbon"][EntityProperty["Element", "AtomicMass"]] + Entity["Element", "Oxygen"][ EntityProperty["Element", "AtomicMass"]]), "AtomicMassUnits"]

In[5]:=\[Omega] = Sqrt[QuantityMagnitude[Quantity[1.86, "Kilonewtons"/"Meters"],"AtomicMassUnit"/"Femtoseconds"^2]/m]

In[6]:=\[HBar] = QuantityMagnitude[Quantity[1., "ReducedPlanckConstant"], "AtomicMassUnit"*"Picometers"^2/"Femtoseconds"]

位移概率密度函数由

给出.

In[7]:=\[Rho][x_, t_] = FullSimplify[ComplexExpand[Conjugate[\[Psi][x, t]] \[Psi][x, t]]]

由于是概率分布，对所有

，ρ 在实数上的积分都为 1.

‍‍ In[8]:=Chop[Integrate[\[Rho][x, t], {x, -\[Infinity], \[Infinity]}]]

可视化概率密度随时间变化的情况.

‍In[9]:=Animate[Plot[\[Rho][x, t], {x, -25, 25}, PlotRange -> {0, .16}, PlotTheme -> "Detailed", FrameLabel -> {Row[{x, RawBoxes@RowBox[{"(", "\"pm\"", ")"}]}," "], None}, LabelStyle -> Larger,PlotLegends ->Placed[{Row[{HoldForm[\[Rho]][x,Quantity[NumberForm[t, {2, 1}], "Femtoseconds"]], RawBoxes@RowBox[{"(", SuperscriptBox["\"pm\"", -1], ")"}]}, " "]}, Above]], {t, 0., 5.7, ImageSize -> Small}, AnimationRate -> 1, SaveDefinitions -> True, Alignment -> Center]

示例2:分析轿车的声学特征模

计算一个 Mini 轿车近似截面的声学特征值和特征函数.

导入一个截面图象.

‍

‍In[1]:=img = Import[ "http://upload.wikimedia.org/wikipedia/commons/d/d3/Mini_cross_\ section.jpg"]

Out[1]:=

使用蒙版工具创建边界图形.

In[2]:=

离散化图形.

In[3]:=bdr = BoundaryDiscretizeGraphics[boundary]

Out[3]:=

计算截面的六个特征值和特征函数.

In[4]:={vals, funs}=NDEigensystem[{-Laplacian[u[x,y], {x, y}]},u[x, y], {x, y}\[Element] bdr,6];

查看特征值.

In[5]:=vals

Out[5]:=

将轿车截面的第二个特征函数可视化.

‍In[6]:=Show[img,ContourPlot[funs[[2]],{x, y} \[Element] bdr, Axes -> None, Frame -> None, AspectRatio -> Automatic, ColorFunction -> Function[f, {Opacity[0.75], ColorData["TemperatureMap"][f]}]], ImageSize -> Automatic]

Out[6]:=

示例3:创建球内的拉普拉斯算子特征函数图集

定义一个三维拉普拉斯算子.

In[1]:=\[ScriptCapitalL] = -Laplacian[u[x, y, z], {x, y, z}];

设定齐次狄利克雷边界条件.

In[2]:=\[ScriptCapitalB] = DirichletCondition[u[x, y, z] == 0, True];

求球内最小的十六个特征值和特征函数.

In[3]:=\[CapitalOmega] = Ball[{0, 0, 0}, 2];

In[4]:={vals, funs} = DEigensystem[{\[ScriptCapitalL], \[ScriptCapitalB]}, u[x, y, z], {x, y, z} \[Element] \[CapitalOmega], 16];

特征值以 BesselJZero 的形式给出.

In[5]:=vals[[1]] // TraditionalForm

生成特征函数图集.

‍In[6]:=Grid[Partition[ParallelTable[DensityPlot3D[ funs[[i]]// N // Evaluate, {x, y, z} \[Element] \[CapitalOmega], Boxed -> False, Axes -> False, ColorFunction -> Hue, Method -> {"ShrinkWrap" -> True}, ImageSize -> 125], {i, 16}], 4]]

Out[6]:=