偏微分を積分しよう
このワークシートはMath by Codeの一部です。
前回は、微分方程式の入門としてのかんたんな変数分離を扱いました。
今回は、方程式のタイプに応じた解法となる、基本的な視点をふやしてみよう。
全微分df=0。(6x+3y+5)dx+(3x-4y +3)dy=0
全微分df=0。(x^3+3xy^2)dx+(3x^2y+y^3)dy=0
１．完全形と積分因子
＜完全形p(x,y)dx+q(x,y)dy=0＞
p(x,y)dx+q(x,y)dy=0の型が、
関数f(x,y)の全微分を使って、
df(x,y)
=∂f/∂x dx +∂f/∂y dy
=p(x,y) dx+q(x,y) dy＝０
と表せる関数ｆがあるとき完全型という。
これは、関数ｆが連続な偏導関数を持ちｆが一定なら全微分df=0となるという意味からきている。
完全形の解は陰関数f(x,y)=C(定数）。完全形が∂p/∂y=∂q/∂xは同値。
∂f/∂x＝p を両辺積分して、f= ∫p dx+c(y)
∂f/∂y＝q を両辺積分して、f= ∫q dy+d(x)
これから、∫pdx+c(y) = ∫qdy+d(x)= f となる、fを定める。
（例）
微分方程式dy/dx=-(6x+3y+5)/(3x-4y +3)
式変形した(6x+3y+5)dx+(3x-4y +3)dy=0をf(x,y)の全微分形式とみなす。
p=6x+3y+5, q=3x-4y +3とおくと、∂p/∂y=∂q/∂x=2だから、完全形。
∫(6x+3y+5) dx +c = 3x2+(3y+5)x+c(y) =3x2+5x+3xy+c(y) =f
∫(3x-4y +3) dy +d = (3x+3)y-2y2+d(x) = d(x) +3xy-2y2 +3y =f
この２式を満たす一般解はf(x,y)=3x2+5x+3xy-2y2 +3y=C
発想は少し高度かもしれないけれど、計算はとても簡単だね。
(別解）
-(6x+3y+5)=3x-4y+3=0の解(x,y) =(-29/33, 1/11)を使って、
x=X-29/33, y=Y+1/11とおいて定数項を消して、変数分離型にもちこむ。
dY/dX=dy/dx=-11(6x+3y+5)/11(3x-4y+3)=-(66(X-29/33)+33(Y+1/11)+55)/(33(X-29/33)-44(Y+1/11)+33)
=-(66X-58+33Y+3+55)/(33X-29-44Y-4+33)=-(66X+33Y)/(33X-44Y)=-(6X+3Y)/(3X-4Y)
u=Y/XとおくとdY/dX=-(6+3(Y/X))/(3-4(Y/X))=(6+3u)/(4u-3)。
Y=Xuの両辺微分dY/dX=(Xu)'=u+X・du/dXからu+X・du/dX=(6+3u)/(4u-3)。
X/dX={(6+6u-4u2)/(4u-3)}/duから(4u-3)/2(3+3u-2u2) du=1/X dXと変数分離。
から、 となり、 だから、
3(x+29/33)2+3(x+29/33)(y-1/11)-2(y-1/11)2=C (C=0の場合も含む）
式を整理すると、
3x2+5x+3xy-2y2 +3y+D=C（D=Cの場合も含む）
発想は単純だけで、計算の手間が増えるね。
（例）
微分方程式dy/dx=-(x3+3xy2)/(3x2y+y3)
式変形した(x3+3xy2)dx+(3x2y+y3)dy=0をf(x,y)の全微分形式とみなす。
p=x3+3xy2, q=3x2y+y3とおくと、∂p/∂y=∂q/∂x=6xyだから、完全形。
∫(x3+3xy2) dx+c= 1/4x4+1/2(3y2)x2+c(y) =1/4(x4+6x2y2+4c(y)) =f
∫(3x2y+y3) dy+d= 1/2(3x2)y2+1/4y4+d(x) =1/4( d(x) +6x2y2+y4) =f
この２式を満たす一般解はf(x,y)=1/4(x4+6x2y2+y4)=C
＜積分因子＞
Pdx+Qdy=0が完全微分形でなくても、適当な式Mをかけて完全型になるときMを積分因子という。Mの求め方としては次のやり方がうまくいくときがある。
N=(Py-Qx)/Qがｘだけの関数なら、M(x)=exp(∫Ndx)
N=(Qx-Py)/Pがｙだけの関数なら、M(y)＝exp(∫Ndy)
他にP,Qのx、yの指数を見比べるMの指数を考えるなど。
（例）
微分方程式dy/dx=-(1-xy)/(xy-x^2)
(1-xy)dx+(xy-x^2)dy=0に適当なMをかけて完全形にするためのMは？
M＝1/xとして、p=1/x-y, q=y-xとおくと、∂p/∂y=∂q/∂x=-1だから、完全形になった。
∫(p) dx+c=log|x|+c(y) =f(x)
∫(q) dy+d= 1/2y2-yx+d(x)= =f(x)
この２式を満たす一般解はf(x,y)=log|x|+1/2y2-xy=C
積分因子でdf=0。(1-xy)dx+(xy-x^2)dy=0
２．実装
質問：全微分形式の微分方程式をコードで解くにはどうしたらよいでしょうか。
sympyは陽関数の形式でのy'=f(x,y)の形で情報を与えると、陽関数の形式で返します。
全微分形式では、x、yのどちらが独立変数ということはないのですが、yがｘの関数としてかいてみます。
y = Function('y')
dy/dxをDerivative(y(x), x)とかきましょう。
[IN]=======================================
from sympy import symbols, Eq, Function, dsolve
x, y = symbols('x y')
y = Function('y') # y=f(x)
def dydxEq(fx):
deq = Eq(y(x).diff(x), fx)
return dsolve(deq, y(x),implicit=True)
dydxEq(-(6*x+3*y(x)+5)/(3*x-4*y(x) +3))
[OUT]=====================================
[Eq(y(x), 3*x/4 - sqrt(C1 + 35937*x**2 + 63162*x)/132 + 3/4),
Eq(y(x), 3*x/4 + sqrt(C1 + 35937*x**2 + 63162*x)/132 + 3/4)]
すごい係数になってしまい、陰関数で返すに指定しても陽関数で返されてしますます。
そこで、全微分形式で、f(x,y)=∫pdx + c(y)= ∫qdy + d(x)の考え方で、f(x,y)=Cを求めてみましょう。
fpdxと∫qsyの共通項はxとyの両方の変数をもつ項です。
だから、集合で考えると
A=∫pdx＝ｘのある項、
B＝∫qdy=ｙのある項、
A∩B=xとyの両方ある項
だから、解の方程式はA∪B-A∩B=Cです。
#[IN]Python===============================
#全微分形式 p(x,y) dx + q(x,y) dy=0の解
from sympy import symbols, Eq, Function, dsolve, integrate
import sympy
def perfectDiff(p_xy,q_xy):
# 積分変数を定義
x, y ,C= symbols('x y C')
#p,qを積分して展開する
px = integrate(p_xy, x).expand()
qy = integrate(q_xy, y).expand()
#多項式を合併する
A=set(list(px.args))
B=set(list(qy.args))
poly=list(A.union(B))
return Eq(sum(poly),C)
x, y ,C= symbols('x y C')
p = 6*x+3*y+5
q = 3*x-4*y+3
perfectDiff(p,q)
#[OUT]==================================
#[IN]========================
p=x**3+3*x*y**2
q=3*x**2*y+y**3
perfectDiff(p,q)
#[OUT]======================
#====================
p=1/x-y
q=y-x
perfectDiff(p,q)
#=====================
＜勾配(grad)＞
惑星Jの東経ｘ、北緯ｙの地点の山の高さをf(x,y)としましょう。
ｘを固定したときの山の断面でのｙ方向の勾配が∂/f∂x
ｙを固定したときの山の断面でのｘ方向の勾配が∂f/∂yです。
山の斜面での最大傾斜線の傾きが全微分dfです。
df = ∂/f∂x dx+ ∂f/∂y dy
x,y地点での微小ベクトルdr=(dx, dy)にたいして、
（∂/f∂x, ∂f/∂y)を勾配grad f演算といいます。
df = grad f drとベクトルの内積の結果
質問：ｆをｚ方向として、geogebraで３Dで山と各点でのdf＝0の確認をするにはどうしたらよいでしょうか。
たとえば、冒頭の全微分の場合は、先に曲面の方程式を与えましょう
f(x,y)=3x2+5x+3xy-2y2 +3y
これが山ｆです。
Derivative(f, x)とすると、
∂/f∂x＝6x+3y+5が計算されます。これをp(x,y)とします。
Derivative(f, y)とすると、
∂f/∂y=3x-4y +3が計算されます。これをq(x,y)とします。
山を水平に切るために、スライダーcを作ります。
z=cの数式だけで、高さｃの水平面rができます。
IntersectionPaths(f,r)を使うと、交線eq1が図でかけます。
これって、山の等高線になっているね。
f(x,y)=cと同等な数式を入れると、山を切る交線としての曲線eq2ができますね。
これは、ｚ座標が表面に出てこないので、ｘｙ平面（グラフィックビュー）の曲線になる。
曲線eq2にあるA点でベクトルdr=(dx,dy)を作り、それを動かすためのｘ座標をスライダーaxで作ります。
Solution(f(ax, y)=c,y)とすると、Aのy座標のリストysが得られます。
１つの交点で十分なので、ay=ys(1)とすると、１つめの交点のy座標を求めてくれます。
そこで、Aを(ax,ay)で定義しましょう。
m(x,y)=-p(x,y)/q(x,y)とおき、C=(ax+1, ay+m(ax,ay))とすることで、
dr=Vector(A,C)で接線の方向ベクトルとしての増分ベクトルdrができるね。
どうようにして、勾配ベクトルgrad fをつくります。
D=(ax+p(ax,ay), ay+q(ax,ay))として、gradf=Vector(A,D)とします。
内積gradf *dr＝０なら、
df = ∂/f∂x dx+ ∂f/∂y dy＝０
となります。
グラフィックビューで２つのベクトルが直交するのが確認できますね。
また、gradf ベクトルが相当でかいので、山というより断崖絶壁のような勾配になっているのことも
感じられるでしょう。
このように、f(x,y)を基点にして、数式を連携させると、
f(x,y)の式を変えるだけで、すべてが連動して変わりますね。