ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

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

この例では、波動方程式の例を使用して Symbolic Math Toolbox™ で汎関数微分を使用する方法を説明します。両端が固定された弦の波動方程式を、汎関数微分を使用して解きます。汎関数微分は、汎関数の、自身が依存する関数についての微分です。Symbolic Math Toolbox™ は関数 functionalDerivativeを使用して汎関数微分を実装します。

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

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

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

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

T=0Ld2(ddtS(x,t))2dx.

位置エネルギーは

V=0Lr2(ddxS(x,t))2dx,

ここで r は張力です。

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

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

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

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

Lr2x2 S(x,t)-Ld2t2 S(x,t)=0

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

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

2t2 S(x,t)v2=2x2 S(x,t)

変数分離法を使用して方程式を解きます。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) = 

2t2 V(t)v2V(t)=2x2 U(x)U(x)

tmp = children(eqn2);

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

syms C
eqn3 = tmp(1) == C
eqn3 = 

2t2 V(t)v2V(t)=C

eqn4 = tmp(2) == C
eqn4 = 

2x2 U(x)U(x)=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) = -2C5sinh(Ctv)
U(x) = simplify(U(x),'Steps',50)
U(x) = 2C8sinh(Cx)

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

p1 = setdiff(symvar(U(x)),sym([C,x]))
p1 = C8
p2 = setdiff(symvar(V(t)),sym([C,v,t]))
p2 = C5

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

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

-k2π2L2

param = k
cond = C801kkZ
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)]);

振動の振幅はパラメーター p1p2 が決定します。単純化のために、p1p21 に設定します。

S = subs(S,[p1 p2],[1 1]);
S = simplify(S,'Steps',50)
S(x, t) = 4sin(πkt)sin(πkx)

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

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)=k=nmAksin(πkt)sin(πkx),

ここで Ak は任意の定数です。

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)=k=nmAksin(πkt)sin(πkx)Ak の値は、初期速度の条件を指定して計算できます。

ut(x,t=0)=Ft(x,0).

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