月別アーカイブ: 2019年4月

プログラミングと数値計算, 物理学

質点と壁との反発を表す運動方程式

コメントする

質点が壁に衝突し、反発することを数式で表現します。

  1. 壁の定義
  2. 運動方程式の導出
  3. プログラム
  4. 実行結果
  5. 参考文献

壁の定義

壁とは、壁の法線方向に平行な質点の速度成分を反転させるデバイスである。
と定義します。

運動方程式の導出

まず議論を簡単にする為に一次元の運動について考え、その後多次元の運動について定式化をしていきます。

1次元の壁との衝突を表す運動方程式

ポテンシャル\(V(x)\)保存力下の質点の運動を考えます。
壁が
\(f(x,t)=0\)
で表現されているとします。
壁との衝突では位置は連続、速度は不連続な振る舞いをすると考えると、壁に衝突した時、力はデルタ関数の振る舞いをしていることが予想できます。
よって、未知の定数\(c\)を用いて、

と書くことが出来ます。
ここで、時刻\(t＝t_i\)は、
\(
f(x(t),t)=0~~\to~~t=t_i
\)
を満たす時刻(壁との衝突時刻)です。\(f(x,t)\)の時間微分は、

と書くことが出来ます。
\(c\)を定めるために式(2)の両辺を、時刻\(t=t_j\)周りを微小時間\(\Delta\)で積分します。すると

を得ます。ここで、保存力の時間変化は連続であることを仮定します。すなわち、

がいかなる時刻\(t\)で成立するとします。計算を進めると、

という結果が導かれます。ここで、\(\delta_{i,j}\)はクロネッカーのデルタを表します。式(6)を式(2)の右辺に代入し、\(c\)を消去すると、運動方程式

を得ます。式(7)はこのままでは解くことが出来ません。なぜなら、未来の時刻の速度\(v(t_i+0)\)が含まれているからです。これをどうにかするには、新たな条件式、すなわち反発前後の条件式が必要になります。

壁と衝突する場合、質点の衝突前後の速度\(v(t_i-0), v(t_i+0)\)と、壁の衝突前後の速度\(v_\text{w}(t_i-0), v_\text{w}(t_i+0)\)の間には、反発係数\(e\)を用いて

の関係があると予想します。今、壁の質量が無限大であり、衝突前後で速度変化がない場合を考えると、壁の速度は、質点との衝突によって変化しないと考えます。すなわち

が成立していると考えます。すると、衝突後の質点の速度は、\(e\)を用いて

と書けます。よって、式(7)に代入して、

を得ます。

最後に、位置\(x(t_i)\)における壁の速度\(v_\text{w}(t_i)\)を考えましょう。衝突の前後のごく短い時間では、質点の動きは壁に追従すると考えます。すなわち、

が成り立っているとします。衝突時には質点の位置\(x\)と衝突位置は同じであることを注記します。式(12)を書き換えると、衝突が起こる時刻\(t=t_i\)の周りで

が成立しています。このことから、衝突前後のごく短い時間では

が成立します。よって、

から、壁の速度

を得ます。よって、壁との衝突を記述する運動方程式は

となります。

多次元の壁との衝突を表す運動方程式

壁

と質点が衝突することを考えます。運動方程式は

であり、1次元の場合と同様に時刻\(t=t_j\)周りの微小時間で積分して

を得ます。式(20)に代入すれば運動方程式

を得ます。次の節で述べる結果(多次元の反発)を先取りすると、式(22)は

と変形することが出来ます。ここで、\(e\)は反発係数、\(\mathbf{v}_\text{w}(t)\)は壁の速度ベクトル、\(\mathbf{n}\)は壁の法線ベクトルであり、\(\mathbf{n}\)は

と表されます。
実際に数値計算を行う際には、衝突時刻\(t_i\)と位置\(\mathbf{r}=\mathbf{r}_i\)を求めた後、

に従って速度ベクトルを変更すると良いと思います（根拠は特にありません）。

多次元の反発


衝突時刻\(t=t_i\), 点\(\mathbf{r}(t_i)\)で質点が\(f(\mathbf{r},t)=0\)で表される壁に衝突することを考えます。
衝突前、後の質点の速度ベクトルをそれぞれ\(\mathbf{v}(t_i-0), \mathbf{v}(t_i+0),\)と置きます。
位置\(\mathbf{r}_i\)における壁の単位法線ベクトルを\(\mathbf{n}, (\mathbf{n}^2=1)\)と書くと、

と書くことが出来ます。ここで、\(c\)は未知の定数でこれから定めていきます。
衝突では壁の法線方向の速度成分のみが変化すると仮定しているので、

が成立します。ここで、\(e\)は反発係数で\(0\leq e\leq 1\)である(\(e=0\):完全非弾性衝突、\(e=1\):完全弾性衝突。
また、\(v_\perp, v_{\text{w}\perp}\)は\(\mathbf{v}, \mathbf{v}_\text{w}\)の、壁の法線方向の成分であり、

と書きあらわすことが出来ます。
式(27)の\(c\)を求めましょう。式(27)の両辺に\(\mathbf{n}^\mathsf{T}\)を掛けて

のように得ます。よって衝突後の質点の速度は

と表すことが出来ます。壁が\(f(\mathbf{r},t)=0\)を満たす線と表現されていれば、時刻\(t=t_i\)で\(\mathbf{n}\)は

と求める事が出来ます。
続いて壁面の速度\(\mathbf{v}_\text{w}(t)\)を考えます。1次元の場合と同様に衝突前後のごく短い時間では質点の動きは壁に追従すると予想します。すなわち、

が成立すると考えます。ここで、式(35)の\(\mathbf{r}(t)\)は質点の位置ベクトルです。すなわち、衝突が起こる時刻\(t=t_i\)の周りで

が成立しており、衝突のごく短い時間では

が得られます。よって

を得ますので、壁の速度

を得ることが出来ます。

プログラム

2次元平面の運動に対するFortran90のプログラムはこちらです。
時間発展は、刻み幅制御陽的ルンゲクッタ法
衝突時刻を求める際の根の探索には、Anderson-Björk’s法を用いています。

壁の関数(サブルーチン fw)とその偏微分(サブルーチン pwf)は手で入力しています。

壁との衝突は、関数の”符号が変わった時”で判定しているので、法線の方向はどちらでも構いません。
すなわち、
\(
f(x,y,t)=\pm g(x,y,t)
\)
は同じです。

プログラム

gnuplot用のスクリプト

実行結果

プログラムを実行した結果です。
適当な壁を定義して実行しています。

重力なしの運動、動かない壁

楕円の焦点から放たれた質点の軌跡

\(
\begin{align}
x(0)&=4,~~x'(0)=(\text{const}) \\
y(0)&=0,~~y'(0)=(\text{const}’)
\end{align}
\)

壁
\(
\begin{gather}
f(x,y,t)=\left(\frac{x}{l_x}\right)^2+\left(\frac{y}{l_y}\right)^2-1=0,\\
l_x=5,~l_y=3
\end{gather}
\)

重力下の運動、動かない壁


\(
\begin{align}
x(0)&=0,~~x'(0)=2 \\
y(0)&=0,~~y'(0)=7
\end{align}
\)

壁
\(
\begin{align}
f(x,y,t)=y – \sin(x)=0
\end{align}
\)

重力下の運動、動く壁


\(
\begin{align}
x(0)&=4,~~x'(0)=0 \\
y(0)&=4,~~y'(0)=5
\end{align}
\)

壁
\(
\begin{gather}
f(x,y,t)=y – \sin(kx+\phi_x)\sin(\omega t+\phi_t)=0,\\
k=1,~w=1,~\phi_x=-\pi/2,~\phi_t=0
\end{gather}
\)

参考文献

3. 壁との衝突 -物理学の見つけ方

4. 動く壁との衝突 -物理学の見つけ方

物理学

量子テレポーテーションのざっくりとした説明

コメントする

量子テレポーテーションとは、量子状態を異なる場所に移動させる方法です。

良くある間違いとして、
・光速を超えて情報通信が出来る(EPRパラドックス)
・物体そのものが移動する(SFのテレポーテーション、ど〇でもドア等)
・物体のコピーが生成される
が挙げられます。上の3点は正しくありません。

本稿の内容が間違えていたら申し訳ありません。

量子テレポーテーションとは

量子テレポーテーションとは、量子もつれを利用して「地点Aの量子状態Xを地点Bに移すこと」です。
量子力学の基本として、未知である任意の量子状態に対し、それと全く同じ複製を作る事は不可能です(量子複製不可能定理 -wikipedia)。
量子状態の複製は出来ませんが、量子状態の移動は行うことが出来ます。これが量子テレポーテーションなのです。
重要なのは、移動するのは物質ではなく、量子状態であることが重要です。

量子テレポーテーションのアイデア

・送りたい側　
　地点Cから送られてきたもつれた量子1’と送りたい量子状態Sを相互作用させる→弱く観測
・受け取り側　
　地点Cから送られてきたもつれた量子2’に対し、弱く観測した結果に基づいた演算をする。

  1. 地点Cで量子1, 2をもつれ(相互作用)させる
    → 量子1,2は区別できないので、量子1′,2’と書くことにする。
  2. 量子1’を地点A、量子2’を地点Bに届ける。
  3. 地点Aで送りたい量子状態Sと量子1’を相互作用させて”弱く”観測する。
    →送りたい量子状態Sはこの時点で量子1’と区別が出来なくなる。
  4. 地点Aの観測結果を古典的な光の速さで地点Bに伝える。
  5. 地点Bの量子2’に、地点Aの観測結果に依存した操作をする。
  6. 地点Bの量子2’は量子状態Sになっている。

という手順です。

地点Aの観測結果を地点Bに送る、という操作が入るので、意味のある情報(量子状態S)を光速を超えて通信する、ということは出来ないわけです。

また、”弱く”観測するの部分では、量子状態を決めることはしません。しっかり観測してしまうと確定的な情報を得てしまい、ダメだそうです。

物体の移動をともなわないにもかかわらず、”テレポーテーション”と呼ばれる由来は、量子が区別できないことに由来します。
区別することができないので、地点Aの量子状態Sと、全てが終わった後の地点Bに存在する量子状態Sは同じものとして捉えられ、あたかもテレポートしたように見える訳です。

同種の粒子が2つあり、それらが同じ量子状態であるのならば、どちらがオリジナルなのか区別する方法はありません。

2017年には中国のグループが、1000km程度の地上-衛星間の量子テレポーテーションに成功した、との報告があります。
Ji-Gang Ren et al, “Ground-to-satellite quantum teleportation”, arXiv:1707.00934 [quant-ph], (2017) https://arxiv.org/abs/1707.00934

参考文献

量子テレポーテーション https://www.m-nomura.com/st/qteleport.html

ーーー観測するって何なんでしょうね。

プログラミングと数値計算, 数学

Hyper-dual numbersによる二階偏微分の計算

コメントする

Hyper-dual numbersと呼ばれる、実数を拡張した考えを導入すると導関数が計算できます。

あらかじめ、引数がHyper-dual numbersである時の数々の関数の定義を実装しておけば、その関数の組み合わせで作られる関数の一階導関数、二階導関数のほぼ厳密な答えを得ることが出来ます。
この考えに従う関数の微分方法は、Forward型の自動微分と呼ばれます。

  1. Dual number
    1. Dual numberと導関数
  2. Hyper-dual numbers
    1. Hyper-dual numbersの演算
    2. Hyper-dual numbersのプログラム
    3. ヘッセ行列
  3. 参考文献

Dual number

Dual numberという考えがあります[2,3,5]。

これは実数を拡張する、という考えで複素数に似た考え方です。

良く知られている実数の拡張方法の一つは、複素数です。
通常の実数に\(i^2=-1\)を満たす\(i\)という数を付加するのが複素空間です。

拡張の方法は何も\(i\)だけではありません。
例えば、\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす\(\epsilon\)という数を付加することもできます。

この\(\epsilon\)は複素数ではありません。
複素平面上でこの性質を満たす数は無いことからも、この\(\epsilon\)を追加するということは新しい方向への数の拡張です。
実際、複素平面上で\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす数があるのかを探しても
\(
\begin{align}
(a+ib)^2=a^2-b^2+i2ab=0
\end{align}
\)
であるので、これを満たすのは\(a=b=0\)しか存在せず、複素平面上の数ではないことが分かります。

そんな新しい方向\(\epsilon\)の平面で定義された数をDual number(二重数、または双対数)と呼びます。

二重数\(a\)は実数部と非実数部から構成されており、
\(
a=a_0+a_1\epsilon
\)
と書くことが出来ます。ここで、\(a_0, a_1\)は実数であり、 \(\epsilon\)は、虚数単位\(i\)に倣って二重数単位とでも名付けておきましょう。

Dual numberと導関数

二重数は関数の導関数と大きな関係がある事を示しましょう。
テーラー展開を行います。
関数\(f\)の\(x\)周りの展開は刻み幅\(\Delta\)を用いて
\(\displaystyle
f(x+\Delta)=f(x)+f'(x) \Delta+\frac{1}{2!}f^{\prime\prime}(x)\Delta^2+\cdots
\)
と記述することが出来ます。もし仮に\(\Delta\)が純非実数だとしましょう。すなわち、\(\Delta=h\epsilon\)とします。ここで、\(h\)は実数です。テーラー展開の式に代入すれば、
\(
\begin{eqnarray}
f(x+h\epsilon)&=& f(x)+f'(x) h\epsilon+\frac{1}{2!}f^{\prime\prime}(x){h}^2\epsilon^2+\cdots \\
&=& f(x)+f'(x) h\epsilon
\end{eqnarray}
\)
となるわけです。注記しますが、上の式は\(h\)の1次で打ち切っているのではなく、厳密にイコールが成り立っています。
この式が言っているのは、関数\(f(a), (a\text{は二重数})\)を計算し、その非実数部を\(h\)で割ると関数の導関数になっているということです。

二重数の非実数部を取り出す関数\(\text{Dual}\)を定義します。
すると、導関数は
\(\displaystyle
f'(x)=\frac{1}{h}\text{Dual}(f(x+h\epsilon))
\)
として得られます。
言葉で書けば、二重数空間に拡張した\(f(x+h\epsilon)\)を計算すると、その非実数部に導関数が現れる、ということです。

二重数のプログラムはJeffrey Fikeさんによる[4]にありますので、そちらをご参考ください。

Hyper-dual numbers

さて、ここまでで二重数の概念を簡単に説明し、新しい実数の拡張を行いました。
Dual Numberのままでは高階導関数は得られません。なぜなら、導関数の二次以降は二重数の性質\(\epsilon^2=0\)によってゼロになるからです。

そこで、若干工夫します。
Dual Numberを二種類用意することで二階微分を得ることが出来ます[2,3]。
すなわち、
\(\begin{gather}
a=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
\epsilon_1^2=\epsilon_2^2=(\epsilon_1\epsilon_2)^2=0,~\epsilon_1\ne 0,~~\epsilon_2\ne 0,~~\epsilon_1\epsilon_2 \ne 0
\end{gather}
\)
と実数を拡張します。この様に拡張した数\(a\)をHyper-dual numbersと呼びます[2,3]。

Hyper-dual numbersの演算

Hyper-dual numbersである
\(
\begin{align}
a&=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
b&=b_0+b_1\epsilon_1+b_2\epsilon_2+b_3\epsilon_1\epsilon_2
\end{align}
\)
を用意します。和、積、商はそれぞれ
\(
\begin{align}
a+b&=(a_0+b_0)+(a_1+b_1)\epsilon_1+(a_2+b_2)\epsilon_2+(a_3+b_3)\epsilon_1\epsilon_2\\
ab&=a_0b_0+(a_0b_1+a_1b_0)\epsilon_1+(a_0b_2+a_2b_0)\epsilon_2+(a_0b_3+a_1b_2+a_2b_1+a_3b_0)\epsilon_1\epsilon_2\\
\frac{1}{a}&=\frac{1}{a_0}-\frac{a_1}{a_0^2}\epsilon_1-\frac{a_2}{a_0^2}\epsilon_2+\left(-\frac{a_3}{a_0^2}+\frac{2a_1a_2}{a_0}\right)\epsilon_1\epsilon_2
\end{align}
\)
と定義されます[2]。任意の関数は
\(
f(x)=f(x_0)+x_1f'(x_0)\epsilon_1+x_2f'(x_0)\epsilon_2
+\left(x_3f'(x_0)+x_1x_2f^{\prime\prime}(x_0)\right)\epsilon_1\epsilon_2
\)
として計算することが出来ます。なので、この結果から\(\epsilon_1\epsilon_2\)の係数として関数の二階微分が得られます。

Hyper-dual numbersのプログラム

実際にプログラムを組みましょう。
Hyper-dual numbersの型を持つ変数はFortran90では定義できません。
なので、構造体を利用して自分で型と、演算を定義します。

Fortran90のプログラムは以下の通りになるかと思います。
基本的な四則演算、基本的な初等関数のHyper-dual numbersの演算をモジュールとして書いています。

上のモジュールと下のメインプログラムを一緒にコンパイルすることにより、関数の二階微分が得られます。

例として
2変数関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)
の偏微分
\(\begin{align}
\frac{\partial f}{\partial x},~~\frac{\partial f}{\partial y},~~\frac{\partial^2 f}{\partial x\partial y}
\end{align}
\)
を得ることを考えます。
\(
f(x+1\epsilon_1,y+1\epsilon_2)
\)
を計算すると
\(
f(x+1\epsilon_1,y+1\epsilon_2)=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2
\)
のように実数係数\(a_0, a_1, a_2, a_3\)が得られます。
すると、
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial y}=a_2 \\
\frac{\partial^2 f}{\partial x\partial y} =a_3
\end{align}
\)
として偏微分が得られます。
ちなみに、二階微分が欲しい場合は
\(
f(x+1\epsilon_1+1\epsilon_2, y)
\)
を考えると
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial x}=a_2 \\
\frac{\partial^2 f}{\partial x^2} =a_3
\end{align}
\)
として得られます。\(a_1,~a_2\)はどちらを採用しても構いません。

プログラムでは変数の型Hyperdualを持つ入力変数を\(\text{xH,yH}\),出力を\(\text{rH}\)と置いています。

program main
  use Hyperdualmod
  implicit none
  type(Hyperdual)::xH,yH,rH
 
  xH%x0 = 0.3d0 ! real part
  xH%x1 = 1d0   ! unreal part \epsilon_1
  xH%x2 = 0d0   ! unreal part \epsilon_2
  xH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  yH%x0 = 0.4d0 ! real part
  yH%x1 = 0d0   ! unreal part \epsilon_1
  yH%x2 = 1d0   ! unreal part \epsilon_2
  yH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  write(6,'(4f23.16)')xH%x0,xH%x1,xH%x2,xH%x3
  write(6,'(4f23.16)')yH%x0,yH%x1,yH%x2,yH%x3  
 
  rH = log(xH*yH**2)*exp(xH)/sqrt(sin(xH)**3+cos(yH)**3)
  !rH = asin(2d0*xH)*acos(yH)/atan(xH*yH)
  !rH = xH**yH
 
  write(6,'(4f23.16)')rH%x0,rH%x1,rH%x2,rH%x3  

  stop
end program main

ヘッセ行列

二階偏微分の計算が出来たので、ヘッセ行列が簡単に計算できます。
ルーチンを作れば、以下の通りになるかと思います。
下のプログラムは3変数関数
\(\displaystyle
f(x,y,z)=\exp(xy)\tan(z)
\)
の\(x=-2,~y=3,~z=1\)におけるヘッセ行列を計算します。

program main
  use Hyperdualmod
  implicit none
  integer::N
  double precision::x,y,z,f
  double precision,allocatable::nabla(:),Hesse(:,:)
  double precision,allocatable::w(:)
  external::func

  N=3

  allocate(nabla(1:N),Hesse(1:N,1:N))
  nabla=0d0
  Hesse=0d0
  allocate(w(1:N))
  w=0d0

  x = -2d0
  y =  3d0
  z =  1d0
  w(1) = x
  w(2) = y
  w(3) = z
  call Hessian(N,w,func,nabla,Hesse)

  f=exp(x*y)*tan(z)
  write(6,'(1e24.16)')nabla(1)
  write(6,'(1e24.16)')y*f
  write(6,'(1e24.16)')nabla(2)
  write(6,'(1e24.16)')x*f
  write(6,'(1e24.16)')nabla(3)
  write(6,'(1e24.16)')exp(x*y)/(cos(z)**2)
  write(6,'(3e24.16)')Hesse(1,1:3)
  write(6,'(3e24.16)')f*y**2, f*(1+x*y), f*y/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(2,1:3)
  write(6,'(3e24.16)')f*(1+x*y), x**2*f, f*x/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(3,1:3)
  write(6,'(3e24.16)')f*y/(cos(z)*sin(z)), f*x/(cos(z)*sin(z)), f*2d0/(cos(z)**2)

  stop
end program main

subroutine func(N,x,f)
  use Hyperdualmod
  implicit none
  integer::N
  type(Hyperdual),intent(in)::x(1:N)
  type(Hyperdual),intent(out)::f
 
  f = exp(x(1)*x(2))*tan(x(3))
 
  return
end subroutine func

subroutine Hessian(N,x,func,nabla,Hesse)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::x(1:N)
  double precision,intent(out)::nabla(1:N)
  double precision,intent(out)::Hesse(1:N,1:N)
  external::func

  integer::i,j,k
  type(Hyperdual),allocatable::xH(:)
  type(Hyperdual)::rH
 
  allocate(xH(1:N))
  do i=1,N
     xH(i)%x0 = 0d0
     xH(i)%x1 = 0d0
     xH(i)%x2 = 0d0
     xH(i)%x3 = 0d0
  enddo
 
  do i=1,N
     xH(i)%x0 = x(i)
     xH(i)%x3 = 0d0
  enddo

  do i=1,N
     do k=1,N
        xH(k)%x1 = 0d0
        xH(k)%x2 = 0d0
     enddo
     xH(i)%x1 = 1d0
     xH(i)%x2 = 1d0
     
     call func(N,xH,rH)
     nabla(i) = rH%x1
     Hesse(i,i) = rH%x3
  enddo

  do i=1,N
     do j=i+1,N
        do k=1,N
           xH(k)%x1 = 0d0
           xH(k)%x2 = 0d0
        enddo
        xH(i)%x1 = 1d0
        xH(i)%x2 = 0d0
        xH(j)%x1 = 0d0
        xH(j)%x2 = 1d0
        call func(N,xH,rH)
        Hesse(i,j) = rH%x3
        Hesse(j,i) = Hesse(i,j)
     enddo
  enddo  
 
  return
end subroutine Hessian

Hyperdual.f90にモジュールを入れ、メインプログラムをmain.f90に入れたとすると、以下の結果を得ます。

> gfortran Hyperdual.f90 main.f90  
> ./a.out
  0.1158128336233602E-01
  0.1158128336233602E-01
 -0.7720855574890680E-02
 -0.7720855574890680E-02
  0.8491012233306163E-02
  0.8491012233306162E-02
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458059E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458058E-01

実行結果の奇数行目はHyper-dual Numberによる計算結果、偶数行目は解析解を表します。
また、6行目までは一階微分、7行目以降はヘッセ行列を表します。

参考文献

[1]関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)
の偏微分\(\frac{\partial^2 f}{\partial x\partial y}\)
の計算
https://www.wolframalpha.com/input/?i=D%5Bln(x*y%5E2)*e%5E(x)%2Fsqrt(sin%5E3(x)%2Bcos%5E3(y)),+y,x%5D

[2]Jeffrey A. Fike and Juan J. Alonso,~”The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations”, 49th AIAA Aerospace Sciences Meeting(2011)
http://adl.stanford.edu/hyperdual/fike_aiaa-2011-886_slides.pdf,
J. A. Fike and J. J. Alonso. The Development of Hyper-Dual Numbers for Exact Second Derivative Calculations. AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting, January 4-7, 2011.http://adl.stanford.edu/hyperdual/Fike_AIAA-2011-886.pdf

[3]JeffreyA.Fike,~”Derivative Calculations Using Hyper-Dual Numbers”, Sandia National Laboratories (2016)
https://www.osti.gov/servlets/purl/1368722

[4]Jeffrey Fike, Aerospace Design Lab, http://adl.stanford.edu/hyperdual/

[5]松本佳彦, 新しい数をつくる, (2018) http://ymatz.net/assets/docs/20180629-jtpa-slide-mod