このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

汎関数微分のチュートリアル

この例では、波動方程式の例を使用して Symbolic Math Toolbox™ で汎関数微分を使用する方法を説明します。この例では、計算をシンボリックに実行して解析結果を取得します。両端が固定された弦の波動方程式を、汎関数微分を使用して解きます。汎関数微分は、汎関数の、自身が依存する関数についての微分です。Symbolic Math Toolbox は関数functionalDerivativeを使用して汎関数微分を実装します。

波動方程式の解は、汎関数微分を応用して求められます。波動方程式とは、弦の振動から電磁波の伝播に至るまでのさまざまな波動を説明する、物理学において重要な方程式の 1 つです。この例で説明している手法は、最速降下問題の求解からシャボン玉の極小曲面の検出まで、変分法のさまざまなアプリに応用することができます。

x = 0 と x = L の 2 点間に吊り下げられた長さ L の弦について考えます。弦は、固有の単位長さ当たりの密度と固有の張力をもちます。後で使用するために、長さ、質量密度、および張力を定数として定義します。簡単にするために、これらの定数を 1 に設定します。

Length = 1;
Density = 1;
Tension = 1;

弦が振動している場合、弦の運動エネルギー密度および位置エネルギー密度は静止状態 S(x,t) からの変位の関数になります。これは位置 x および時間 t によって変化します。d を単位長さ当たりの質量密度とすると、運動エネルギー密度は

$$T=\frac{d}{2}{(\frac{d}{dt}S(x,t))}^2.$$

位置エネルギー密度は

$$V=\frac{r}{2}{(\frac{d}{dx}S(x,t))}^2,$$

ここで、r は張力です。

MATLAB® でこれらの方程式を入力します。長さは正の値なので、この仮定を設定します。この仮定により、simplify は得られる方程式を所定の形式に単純化できます。

syms S(x,t) d r v L
assume(L>0)
T(x,t) = d/2*diff(S,t)^2;
V(x,t) = r/2*diff(S,x)^2;

作用密度 A は T-V です。最小作用の原理によれば、作用は常に最小になります。functionalDerivative を使用して A の汎関数微分を S について求め、それを 0 に等しくすることによって、最小作用の条件を決定します。

A = T-V;
eqn = functionalDerivative(A,S) == 0

eqn(x, t) = 
$$r\hspace{0.17em}\frac{{\partial }^2}{\partial x^2}\mathrm{ }S\left(x,t\right)-d\hspace{0.17em}\frac{{\partial }^2}{\partial t^2}\mathrm{ }S\left(x,t\right)=0$$

simplify を使用して方程式を単純化します。r/d を波の速度 v の 2 乗で置き換えて、方程式を所定の形式に変換します。

eqn = simplify(eqn)/r;
eqn = subs(eqn,r/d,v^2)

eqn(x, t) = 
$$\frac{\frac{{\partial }^2}{\partial t^2}\mathrm{ }S\left(x,t\right)}{v^2}=\frac{{\partial }^2}{\partial x^2}\mathrm{ }S\left(x,t\right)$$

変数分離法を使用して方程式を解きます。S(x,t) = U(x)*V(t) を設定して、位置 x および時間 t への依存関係を分離します。children を使用して、得られた方程式を両辺に分離します。

syms U(x) V(t)
eqn2 = subs(eqn,S(x,t),U(x)*V(t));
eqn2 = eqn2/(U(x)*V(t))

eqn2(x, t) = 
$$\frac{\frac{{\partial }^2}{\partial t^2}\mathrm{ }V\left(t\right)}{v^2\hspace{0.17em}V\left(t\right)}=\frac{\frac{{\partial }^2}{\partial x^2}\mathrm{ }U\left(x\right)}{U\left(x\right)}$$

tmp = children(eqn2);

方程式の両辺はそれぞれ異なる変数に依存していますが、等しくなっています。これが成り立つのは、両辺が定数の場合のみです。両辺が任意定数 C と等しいとして、2 つの微分方程式を得ます。

syms C
eqn3 = tmp(1) == C

eqn3 = 
$$\frac{\frac{{\partial }^2}{\partial t^2}\mathrm{ }V\left(t\right)}{v^2\hspace{0.17em}V\left(t\right)}=C$$

eqn4 = tmp(2) == C

eqn4 = 
$$\frac{\frac{{\partial }^2}{\partial x^2}\mathrm{ }U\left(x\right)}{U\left(x\right)}=C$$

dsolve を使用して、x = 0 および t = 0 における変位が 0 という条件で微分方程式を解きます。Steps オプションを 50 に設定し、simplify を使用して方程式を所定の形式に単純化します。

V(t) = dsolve(eqn3,V(0)==0,t);
U(x) = dsolve(eqn4,U(0)==0,x);
V(t) = simplify(V(t),'Steps',50)

V(t) = $$-2\hspace{0.17em}{C}_1\hspace{0.17em}\sinh \left(\sqrt{C}\hspace{0.17em}t\hspace{0.17em}v\right)$$

U(x) = simplify(U(x),'Steps',50)

U(x) = $$2\hspace{0.17em}{C}_1\hspace{0.17em}\sinh \left(\sqrt{C}\hspace{0.17em}x\right)$$

方程式の定数を求めます。

p1 = setdiff(symvar(U(x)),sym([C,x]))

p1 = $$C_1$$

p2 = setdiff(symvar(V(t)),sym([C,v,t]))

p2 = $$C_1$$

弦は位置 x = 0 と x = L に固定されています。条件 U(0) = 0 が既に存在します。境界条件 U(L) = 0 を適用して C を求めます。

eqn_bc = U(L) == 0;
[solC,param,cond] = solve(eqn_bc,C,'ReturnConditions',true)

solC = 
$$-\frac{k^2\hspace{0.17em}{\pi }^2}{L^2}$$

param = $$k$$

cond = $$k\in \mathbb{Z}\wedge C_1\ne 0\wedge 1\le k$$

assume(cond)

解 S(x,t) は U(x) と V(t) の積です。解を求め、弦の特性値を解に代入して解の最終形式を得ます。

S(x,t) = U(x)*V(t);
S = subs(S,C,solC);
S = subs(S,[L v],[Length sqrt(Tension/Density)]);

振動の振幅はパラメーター p1 と p2 が決定します。単純化のために、p1 と p2 を 1 に設定します。

S = subs(S,[p1 p2],[1 1]);
S = simplify(S,'Steps',50)

S(x, t) = $$4\hspace{0.17em}\sin \left(\pi \hspace{0.17em}k\hspace{0.17em}t\right)\hspace{0.17em}\sin \left(\pi \hspace{0.17em}k\hspace{0.17em}x\right)$$

弦の振動状態は、k のさまざまな値に応じて異なります。任意の時間値 t における最初の 4 つの状態をプロットします。solve によって返される引数 param を使用してパラメーター k を参照します。

Splot(x) = S(x,0.3);
figure(1)
hold on
grid on
ymin = double(coeffs(Splot));
for i = 1:4
    yplot = subs(Splot,param,i);
    fplot(yplot,[0 Length])
end
ylim([-ymin ymin])
legend('k = 1','k = 2','k = 3','k = 4','Location','best')
xlabel('Position (x)')
ylabel('Displacement (S)')
title('Modes of a string')

波動方程式は線形です。つまり、可能な状態のどのような線形結合も波動方程式の有効な解になります。したがって、与えられた境界条件と初期値をもつ波動方程式の完全解は、可能な状態の総和となります

$$F(x,t)=\sum _{k=n}^m{A}_k\sin (\pi kt)\sin (\pi kx),$$

ここで $$A_k$$ は任意の定数です。

symsum を使用して、弦の最初の 5 つの状態の合計を求めます。比較のために、新しい図に前の波形と同じ時刻に得られた波形を表示します。

figure(2)
S5(x) = 1/5*symsum(S,param,1,5);
fplot(subs(S5,t,0.3),[0 Length])
ylim([-ymin ymin])
grid on
xlabel('Position (x)')
ylabel('Displacement (S)')
title('Summation of first 5 modes')

図では、状態の加算によって定性的に異なる波形がモデル化されることがわかります。ここではすべての $$x$$ について初期条件を $$S(x,t=0)=0$$ と指定しています。

方程式 $$F(x,t)=\sum _{k=n}^m{A}_k\sin (\pi kt)\sin (\pi kx)$$ の $$A_k$$ の値は、初期速度の条件を指定して計算できます。

$$u_t(x,t=0)=F_t(x,0).$$

どのような波形も状態の適切な総和により表すことができます。これは、フーリエ級数を使用して弦の振動を表す場合も同様です。