{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4c0175fe-1aa9-47ff-a15d-09a524639525",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 一般超幾何関数と超幾何微分方程式\n",
    "\n",
    "一般超幾何関数が満たす微分方程式も超幾何微分方程式と呼びます。最も簡単な一般超幾何関数\n",
    "$$F(x)=_{3}F_2\\left(a,b,c;d,e;x\\right)$$\n",
    "の時、\n",
    "$$x^2\\,(1-x)\\frac{d^3\\,F(x)}{d\\,x^2}+x\\,((d+e+1)-(a+b+c+3)\\,x)\\,\\frac{d^2\\,F(x)}{d\\,x^2}+(d\\,e-(a+b+c+ab+bc+ca+1)\\,z)\\,\\frac{d\\,F(x)}{d\\,x}-a\\,b\\,c\\,F(x)=0$$\n",
    "が成り立ちます。これを証明していきます。\n",
    "\n",
    "証明の方針としてはガウスの超幾何関数の場合とほぼ同様です。$A_n=\\frac{(a)_n\\,(b)_n\\,(c)_n}{(d)_n\\,(e)_n}$として、\n",
    "$$F(x)=\\sum_{n=0}^{\\infty}\\frac{A_n\\,x^n}{n!}$$\n",
    "から$\\frac{d}{dx}F(x), \\frac{d^2}{dx^2}F(x), \\frac{d^3}{dx^3}F(x)$の級数展開を、係数に$A_{n+1}$が現れる場合と$A_n$が現れる場合と２通り計算します。\n",
    "\n",
    "次に上記の$A_n$の定義から自明に導かれる漸化式$A_{n+1}=\\frac{(a+n)\\,(b+n)\\,(c+n)}{(d+n)\\,(e+n)}\\,A_n$から$(d+n)\\,(e+n)\\,A_{n+1}=(a+n)\\,(b+n)\\,(c+n)\\,A_n$の両辺を展開して、さらに両辺の総和を取ります。その各項を先ほど求めた微分に置き換えていくことでこの微分方程式を求めます。\n",
    "\n",
    "前の記事で紹介したガウスの超幾何関数の超幾何微分方程式を求める過程では上記最後のステップの微分に置き換えるところを、「じっと睨んで」適当な係数をかけることで置き換えがうまくようにしていました。今回はその部分も機械的な計算で求められることをMaximaで示します。これによって、より複雑な一般超幾何関数でも同じ手法で対応する超幾何微分方程式が求められることが分かります。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b1c3235-73a8-4624-91fd-ee08b115cc6f",
   "metadata": {},
   "source": [
    "次のセルから定義や計算が始まりますが、しばらくはガウスの超幾何関数と同じ進行ですから、特に解説はありません。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e19e4961-4f1d-41b4-821e-cb58460e4773",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{0}$}{\\it out1stTerm}\\left({\\it sumexp}\\right):=\\mathbf{block}\\;\\left(\\left[ {\\it ar}:{\\it args}\\left({\\it sumexp}\\right) , \\exp  , {\\it varname} , {\\it initvalue} , {\\it finvalue} \\right]  , \\left[ \\exp  , {\\it varname} , {\\it initvalue} , {\\it finvalue} \\right] :{\\it ar} , {\\it subst}\\left(0 , {\\it varname} , \\exp \\right)+{\\it apply}\\left({\\it sum} , \\left[ {\\it subst}\\left({\\it varname}+1 , {\\it varname} , \\exp \\right) , {\\it varname} , {\\it initvalue} , {\\it finvalue} \\right] \\right)\\right)\\]"
      ],
      "text/plain": [
       "(%o0) out1stTerm(sumexp) := block([ar : args(sumexp), exp, varname, initvalue, \n",
       "finvalue], [exp, varname, initvalue, finvalue] : ar, \n",
       "subst(0, varname, exp) + apply(sum, [subst(varname + 1, varname, exp), \n",
       "varname, initvalue, finvalue]))"
      ],
      "text/x-maxima": [
       "out1stTerm(sumexp):=block([ar:args(sumexp),exp,varname,initvalue,finvalue],\n",
       "           [exp,varname,initvalue,finvalue]:ar,\n",
       "           subst(0,varname,exp)+apply(sum,\n",
       "                                      [subst(varname+1,varname,exp),varname,\n",
       "                                       initvalue,finvalue]))"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out1stTerm(sumexp):=block([ar:args(sumexp),exp,varname,initvalue,finvalue],\n",
    "  [exp,varname,initvalue,finvalue]:ar,\n",
    "  subst(0,varname,exp)+apply(sum,[subst(varname+1,varname,exp),varname,initvalue,finvalue]));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "30db0b79-40c7-4d9b-9608-3c152513b789",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{1}$}F\\left(x\\right)=\\sum_{n=0}^{\\infty }{\\frac{A_{n}\\,x^{n}}{n!}}\\]"
      ],
      "text/plain": [
       "                                     inf       n\n",
       "                                     ====  A  x\n",
       "                                     \\      n\n",
       "(%o1)                         F(x) =  >    -----\n",
       "                                     /      n!\n",
       "                                     ====\n",
       "                                     n = 0"
      ],
      "text/x-maxima": [
       "F(x) = 'sum((A[n]*x^n)/n!,n,0,inf)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F0:F(x)=sum(A[n]/n!*x^n,n,0,inf);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "15298aec-7ece-4014-834a-4080a7c29d9d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{2}$}\\frac{d}{d\\,x}\\,F\\left(x\\right)=\\sum_{n=0}^{\\infty }{\\frac{n\\,A_{n}\\,x^{n-1}}{n!}}\\]"
      ],
      "text/plain": [
       "                                     inf         n - 1\n",
       "                                     ====  n A  x\n",
       "                         d           \\        n\n",
       "(%o2)                    -- (F(x)) =  >    -----------\n",
       "                         dx          /         n!\n",
       "                                     ====\n",
       "                                     n = 0"
      ],
      "text/x-maxima": [
       "'diff(F(x),x,1) = 'sum((n*A[n]*x^(n-1))/n!,n,0,inf)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FD1_1:'diff(F(x),x)=diff(rhs(F0),x);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4afd68c3-40de-4342-b7a8-5a5ef7110463",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{3}$}\\frac{d}{d\\,x}\\,F\\left(x\\right)=\\sum_{n=0}^{\\infty }{\\frac{A_{n+1}\\,x^{n}}{n!}}\\]"
      ],
      "text/plain": [
       "                                      inf           n\n",
       "                                      ====  A      x\n",
       "                          d           \\      n + 1\n",
       "(%o3)                     -- (F(x)) =  >    ---------\n",
       "                          dx          /        n!\n",
       "                                      ====\n",
       "                                      n = 0"
      ],
      "text/x-maxima": [
       "'diff(F(x),x,1) = 'sum((A[n+1]*x^n)/n!,n,0,inf)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SB-KERNEL:REDEFINITION-WITH-DEFUN: redefining MAXIMA::SIMP-UNIT-STEP in DEFUN\n",
      "SB-KERNEL:REDEFINITION-WITH-DEFUN: redefining MAXIMA::SIMP-POCHHAMMER in DEFUN\n"
     ]
    }
   ],
   "source": [
    "FD1_2:lhs(FD1_1)=ev(out1stTerm(rhs(FD1_1)),factorial_expand:true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9f276802-360e-4c57-a650-447a46d158c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{4}$}\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)=\\sum_{n=0}^{\\infty }{\\frac{\\left(n-1\\right)\\,n\\,A_{n}\\,x^{n-2}}{n!}}\\]"
      ],
      "text/plain": [
       "                                 inf                 n - 2\n",
       "                     2           ====  (n - 1) n A  x\n",
       "                    d            \\                n\n",
       "(%o4)               --- (F(x)) =  >    -------------------\n",
       "                      2          /             n!\n",
       "                    dx           ====\n",
       "                                 n = 0"
      ],
      "text/x-maxima": [
       "'diff(F(x),x,2) = 'sum(((n-1)*n*A[n]*x^(n-2))/n!,n,0,inf)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FD2_1:'diff(F(x),x,2)=diff(rhs(FD1_1),x);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8b7b1736-d750-43fd-b153-11d073725b2d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{5}$}\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)=\\sum_{n=0}^{\\infty }{\\frac{n\\,A_{n+1}\\,x^{n-1}}{n!}}\\]"
      ],
      "text/plain": [
       "                                   inf             n - 1\n",
       "                       2           ====  n A      x\n",
       "                      d            \\        n + 1\n",
       "(%o5)                 --- (F(x)) =  >    ---------------\n",
       "                        2          /           n!\n",
       "                      dx           ====\n",
       "                                   n = 0"
      ],
      "text/x-maxima": [
       "'diff(F(x),x,2) = 'sum((n*A[n+1]*x^(n-1))/n!,n,0,inf)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FD2_2:lhs(FD2_1)=ev(out1stTerm(rhs(FD2_1)),factorial_expand:true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3564a51a-2c41-4b25-a5fa-66305273e81a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{6}$}\\frac{d^3}{d\\,x^3}\\,F\\left(x\\right)=\\sum_{n=0}^{\\infty }{\\frac{\\left(n-2\\right)\\,\\left(n-1\\right)\\,n\\,A_{n}\\,x^{n-3}}{n!}}\\]"
      ],
      "text/plain": [
       "                             inf                         n - 3\n",
       "                 3           ====  (n - 2) (n - 1) n A  x\n",
       "                d            \\                        n\n",
       "(%o6)           --- (F(x)) =  >    ---------------------------\n",
       "                  3          /                 n!\n",
       "                dx           ====\n",
       "                             n = 0"
      ],
      "text/x-maxima": [
       "'diff(F(x),x,3) = 'sum(((n-2)*(n-1)*n*A[n]*x^(n-3))/n!,n,0,inf)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FD3_1:'diff(F(x),x,3)=diff(rhs(FD2_1),x);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d66f4c9a-e9b9-418e-92b2-e691c2237021",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{7}$}\\frac{d^3}{d\\,x^3}\\,F\\left(x\\right)=\\sum_{n=0}^{\\infty }{\\frac{\\left(n-1\\right)\\,n\\,A_{n+1}\\,x^{n-2}}{n!}}\\]"
      ],
      "text/plain": [
       "                               inf                     n - 2\n",
       "                   3           ====  (n - 1) n A      x\n",
       "                  d            \\                n + 1\n",
       "(%o7)             --- (F(x)) =  >    -----------------------\n",
       "                    3          /               n!\n",
       "                  dx           ====\n",
       "                               n = 0"
      ],
      "text/x-maxima": [
       "'diff(F(x),x,3) = 'sum(((n-1)*n*A[n+1]*x^(n-2))/n!,n,0,inf)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FD3_2:lhs(FD3_1)=ev(out1stTerm(rhs(FD3_1)),factorial_expand:true);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5650f4a1-58bd-46d5-9577-54ccf48a6c03",
   "metadata": {},
   "source": [
    "$A_n$の定義が${}_3 F_2$では異なるため、漸化式も以下のように少し変わります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "407e29be-910f-426c-824f-87da3ad6aa0e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{8}$}A_{n+1}=\\frac{\\left(n+a\\right)\\,\\left(n+b\\right)\\,\\left(n+c\\right)\\,A_{n}}{\\left(n+d\\right)\\,\\left(n+e\\right)}\\]"
      ],
      "text/plain": [
       "                               (n + a) (n + b) (n + c) A\n",
       "                                                        n\n",
       "(%o8)                 A      = --------------------------\n",
       "                       n + 1        (n + d) (n + e)"
      ],
      "text/x-maxima": [
       "A[n+1] = ((n+a)*(n+b)*(n+c)*A[n])/((n+d)*(n+e))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Acond:A[n+1]=(a+n)*(b+n)*(c+n)/((n+d)*(n+e))*A[n];"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcd3cee3-e8a8-45ea-a795-e5bff744359f",
   "metadata": {},
   "source": [
    "分母を払って両辺を展開します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e15432bd-71d1-4b12-b67f-f5bdb794f0df",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{9}$}n^2\\,A_{n+1}+e\\,n\\,A_{n+1}+d\\,n\\,A_{n+1}+d\\,e\\,A_{n+1}=n^3\\,A_{n}+c\\,n^2\\,A_{n}+b\\,n^2\\,A_{n}+a\\,n^2\\,A_{n}+b\\,c\\,n\\,A_{n}+a\\,c\\,n\\,A_{n}+a\\,b\\,n\\,A_{n}+a\\,b\\,c\\,A_{n}\\]"
      ],
      "text/plain": [
       "       2\n",
       "(%o9) n  A      + e n A      + d n A      + d e A      = \n",
       "          n + 1        n + 1        n + 1        n + 1\n",
       " 3         2         2         2\n",
       "n  A  + c n  A  + b n  A  + a n  A  + b c n A  + a c n A  + a b n A  + a b c A\n",
       "    n         n         n         n          n          n          n          n"
      ],
      "text/x-maxima": [
       "n^2*A[n+1]+e*n*A[n+1]+d*n*A[n+1]+d*e*A[n+1]\n",
       "  = n^3*A[n]+c*n^2*A[n]+b*n^2*A[n]+a*n^2*A[n]+b*c*n*A[n]+a*c*n*A[n]+a*b*n*A[n]\n",
       "            +a*b*c*A[n]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F4:expand(Acond*(n+d)*(n+e));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a08baf22-42af-47ca-904c-54d67ee1e095",
   "metadata": {},
   "source": [
    "両辺に$\\frac{x^n}{n!}$をかけてから$n$について総和を取ります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a8bbcc0e-2c91-4d6b-bdbc-79633efc31bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{10}$}\\sum_{n=0}^{k}{\\frac{\\left(n^2+\\left(e+d\\right)\\,n+d\\,e\\right)\\,A_{n+1}\\,x^{n}}{n!}}=\\sum_{n=0}^{k}{\\frac{\\left(n^3+\\left(c+b+a\\right)\\,n^2+\\left(\\left(b+a\\right)\\,c+a\\,b\\right)\\,n+a\\,b\\,c\\right)\\,A_{n}\\,x^{n}}{n!}}\\]"
      ],
      "text/plain": [
       "        k      2                            n\n",
       "       ====  (n  + (e + d) n + d e) A      x\n",
       "       \\                             n + 1\n",
       "(%o10)  >    -------------------------------- = \n",
       "       /                    n!\n",
       "       ====\n",
       "       n = 0\n",
       "                 k      3                2                                    n\n",
       "                ====  (n  + (c + b + a) n  + ((b + a) c + a b) n + a b c) A  x\n",
       "                \\                                                          n\n",
       "                 >    ---------------------------------------------------------\n",
       "                /                                n!\n",
       "                ====\n",
       "                n = 0"
      ],
      "text/x-maxima": [
       "'sum(((n^2+(e+d)*n+d*e)*A[n+1]*x^n)/n!,n,0,k)\n",
       "  = 'sum(((n^3+(c+b+a)*n^2+((b+a)*c+a*b)*n+a*b*c)*A[n]*x^n)/n!,n,0,k)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F5:factor(sum(F4*x^n/n!,n,0,k));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "868419dd-a7f5-470a-9b4b-6e25e276c85b",
   "metadata": {},
   "source": [
    "この式の右辺や左辺を微分の線型結合で表すことができるか、という問題を考えます。$F(x)$の微分の無限級数の係数を見ると、$A_n$や$A_{n+1}$に$n$の下降階乗冪$n$, $n\\,(n-1)$, $n\\,(n-1)\\,(n-2)$などが掛かった形をしています。これらは下降冪記号$(n)^-_k=n\\,(n-1)\\cdots (n-k+1)$を使うと$(n)^-_1, (n)^-_2, (n)^-_3$と書くことができます。実は一般にはこれらに$(n)^-_0=1$を加えた多項式集合$\\{(n)^-_i | i=0,1,\\dots,k \\}$は$n$の$k$次多項式の基底を成すことが知られています。\n",
    "\n",
    "例えば上記の式の左辺の$\\frac{A_{n+1}\\,x^n}{n!}$の係数は$n$の2次式なので、$\\{(n)^-_i | i=0,1,2 \\}$の要素の1次結合で書くことがきます。その時最高次数の係数は1で確定しますが、1次、0次の係数は未定なので$p, q$と置くと次の式が成り立ちます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f4ba5cdd-21f2-4661-99af-3a38da52cf5c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{11}$}n^2+\\left(e+d\\right)\\,n+d\\,e=q+n\\,p+\\left(n-1\\right)\\,n\\]"
      ],
      "text/plain": [
       "                   2\n",
       "(%o11)            n  + (e + d) n + d e = q + n p + (n - 1) n"
      ],
      "text/x-maxima": [
       "n^2+(e+d)*n+d*e = q+n*p+(n-1)*n"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F6:n^2+(e+d)*n+d*e=(n-1)*n+p*n+q;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "5062f5c5-dae9-469d-bcc7-021812cf8164",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{12}$}n^2+e\\,n+d\\,n+d\\,e=q+n\\,p+n^2-n\\]"
      ],
      "text/plain": [
       "                     2                                2\n",
       "(%o12)              n  + e n + d n + d e = q + n p + n  - n"
      ],
      "text/x-maxima": [
       "n^2+e*n+d*n+d*e = q+n*p+n^2-n"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F7:F6,expand;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5ee6c75e-d82e-4412-a477-8b22c7c146e0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{13}$}-q-n\\,p+e\\,n+d\\,n+n+d\\,e=0\\]"
      ],
      "text/plain": [
       "(%o13)               (- q) - n p + e n + d n + n + d e = 0"
      ],
      "text/x-maxima": [
       "(-q)-n*p+e*n+d*n+n+d*e = 0"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F8:lhs(F7)-rhs(F7)=0;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dafb7592-2492-47c0-9b05-b3e4e11a82de",
   "metadata": {},
   "source": [
    "上の式が$n\\ge 0$の任意の$n$について成り立つので、少なくとも$n=1, n=0$でも成り立つとして$p,q$を求めます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "29eebd4f-85b4-4b08-9d44-24495bb2f7e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{14}$}\\left[ \\left[ p=e+d+1 , q=d\\,e \\right]  \\right] \\]"
      ],
      "text/plain": [
       "(%o14)                    [[p = e + d + 1, q = d e]]"
      ],
      "text/x-maxima": [
       "[[p = e+d+1,q = d*e]]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SOL:solve([ev(F8,n=0),ev(F8,n=1)],[p,q]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07d2b27a-298d-41a0-9d6e-5c9f9b65c3b2",
   "metadata": {},
   "source": [
    "念のためですが、これらの$p,q$をF6式に代入すると右辺と左辺が等しくなるので、F6は恒等式であることがわかりました。これらの$p,q$を使ってF5式の左辺を微分で置換します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "21d61576-fe12-4e28-97ed-95ab9608bbc2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{15}$}x^2\\,\\left(\\frac{d^3}{d\\,x^3}\\,F\\left(x\\right)\\right)+e\\,x\\,\\left(\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)\\right)+d\\,x\\,\\left(\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)\\right)+x\\,\\left(\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)\\right)+d\\,e\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)=\\sum_{n=0}^{\\infty }{\\frac{\\left(n^2+\\left(e+d\\right)\\,n+d\\,e\\right)\\,A_{n+1}\\,x^{n}}{n!}}\\]"
      ],
      "text/plain": [
       "            3                  2                  2                2\n",
       "        2  d                  d                  d                d\n",
       "(%o15) x  (--- (F(x))) + e x (--- (F(x))) + d x (--- (F(x))) + x (--- (F(x)))\n",
       "             3                  2                  2                2\n",
       "           dx                 dx                 dx               dx\n",
       "                                         inf     2                            n\n",
       "                                         ====  (n  + (e + d) n + d e) A      x\n",
       "                            d            \\                             n + 1\n",
       "                     + d e (-- (F(x))) =  >    --------------------------------\n",
       "                            dx           /                    n!\n",
       "                                         ====\n",
       "                                         n = 0"
      ],
      "text/x-maxima": [
       "x^2*'diff(F(x),x,3)+e*x*'diff(F(x),x,2)+d*x*'diff(F(x),x,2)+x*'diff(F(x),x,2)\n",
       "                   +d*e*'diff(F(x),x,1)\n",
       "  = 'sum(((n^2+(e+d)*n+d*e)*A[n+1]*x^n)/n!,n,0,inf)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F9:factor(expand(sumcontract(intosum(FD3_2*x^2+p*FD2_2*x+q*FD1_2)))),SOL;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a45ef5e8-3edc-49bc-bd31-5f89dc1d3cb2",
   "metadata": {},
   "source": [
    "F5式の右辺についても同様の計算を行います。F5式の右辺の$\\frac{A_n\\,x^n}{n!}$の係数が$n$の3次式であることからこの係数式は$\\{(n)^-_i | i=0,1,2,3 \\}$の要素の1次結合で書くことがきます。最高次の係数は1で決まりますが、2, 1, 0次の係数を$p, q, r$としてこれらを求めます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "572554e2-dae0-4c3f-a451-b6028f23897a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{16}$}n^3+\\left(c+b+a\\right)\\,n^2+\\left(\\left(b+a\\right)\\,c+a\\,b\\right)\\,n+a\\,b\\,c\\]"
      ],
      "text/plain": [
       "                3                2\n",
       "(%o16)         n  + (c + b + a) n  + ((b + a) c + a b) n + a b c"
      ],
      "text/x-maxima": [
       "n^3+(c+b+a)*n^2+((b+a)*c+a*b)*n+a*b*c"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F10:part(F5,2,1,1,1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "c1fe3e1e-af55-4c16-8f3b-250ebfce9512",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{17}$}n^3+\\left(c+b+a\\right)\\,n^2+\\left(\\left(b+a\\right)\\,c+a\\,b\\right)\\,n+a\\,b\\,c=r+n\\,q+\\left(n-1\\right)\\,n\\,p+\\left(n-2\\right)\\,\\left(n-1\\right)\\,n\\]"
      ],
      "text/plain": [
       "        3                2\n",
       "(%o17) n  + (c + b + a) n  + ((b + a) c + a b) n + a b c = \n",
       "                                      r + n q + (n - 1) n p + (n - 2) (n - 1) n"
      ],
      "text/x-maxima": [
       "n^3+(c+b+a)*n^2+((b+a)*c+a*b)*n+a*b*c = r+n*q+(n-1)*n*p+(n-2)*(n-1)*n"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F11:F10=n*(n-1)*(n-2)+p*n*(n-1)+q*n+r;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "107a72ac-90ae-47dd-8a2a-c84a238d0be9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{18}$}n^3+c\\,n^2+b\\,n^2+a\\,n^2+b\\,c\\,n+a\\,c\\,n+a\\,b\\,n+a\\,b\\,c=r+n\\,q+n^2\\,p-n\\,p+n^3-3\\,n^2+2\\,n\\]"
      ],
      "text/plain": [
       "        3      2      2      2\n",
       "(%o18) n  + c n  + b n  + a n  + b c n + a c n + a b n + a b c = \n",
       "                                                    2            3      2\n",
       "                                         r + n q + n  p - n p + n  - 3 n  + 2 n"
      ],
      "text/x-maxima": [
       "n^3+c*n^2+b*n^2+a*n^2+b*c*n+a*c*n+a*b*n+a*b*c = r+n*q+n^2*p-n*p+n^3-3*n^2+2*n"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F12:F11,expand;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e59e13dc-d145-4945-87cd-ce7efc8497df",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{19}$}-r-n\\,q-n^2\\,p+n\\,p+c\\,n^2+b\\,n^2+a\\,n^2+3\\,n^2+b\\,c\\,n+a\\,c\\,n+a\\,b\\,n-2\\,n+a\\,b\\,c=0\\]"
      ],
      "text/plain": [
       "                      2              2      2      2      2\n",
       "(%o19) (- r) - n q - n  p + n p + c n  + b n  + a n  + 3 n  + b c n + a c n\n",
       "                                                      + a b n - 2 n + a b c = 0"
      ],
      "text/x-maxima": [
       "(-r)-n*q-n^2*p+n*p+c*n^2+b*n^2+a*n^2+3*n^2+b*c*n+a*c*n+a*b*n-2*n+a*b*c = 0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F13:lhs(F12)-rhs(F12)=0;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c6efd3f2-146d-4ec6-b114-bcfa5e6cd67e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{20}$}\\left[ \\left[ p=c+b+a+3 , q=\\left(b+a+1\\right)\\,c+\\left(a+1\\right)\\,b+a+1 , r=a\\,b\\,c \\right]  \\right] \\]"
      ],
      "text/plain": [
       "(%o20) [[p = c + b + a + 3, q = (b + a + 1) c + (a + 1) b + a + 1, r = a b c]]"
      ],
      "text/x-maxima": [
       "[[p = c+b+a+3,q = (b+a+1)*c+(a+1)*b+a+1,r = a*b*c]]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SOL2:solve([ev(F13,n=0),ev(F13,n=1),ev(F13,n=-1)],[p,q,r]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "f1d050bf-7f4b-45fc-b54d-2c94e9973c05",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{21}$}x^3\\,\\left(\\frac{d^3}{d\\,x^3}\\,F\\left(x\\right)\\right)+c\\,x^2\\,\\left(\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)\\right)+b\\,x^2\\,\\left(\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)\\right)+a\\,x^2\\,\\left(\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)\\right)+3\\,x^2\\,\\left(\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)\\right)+b\\,c\\,x\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)+a\\,c\\,x\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)+c\\,x\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)+a\\,b\\,x\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)+b\\,x\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)+a\\,x\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)+x\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)+a\\,b\\,c\\,F\\left(x\\right)=\\sum_{n=0}^{\\infty }{\\frac{\\left(n^3+\\left(c+b+a\\right)\\,n^2+\\left(\\left(b+a\\right)\\,c+a\\,b\\right)\\,n+a\\,b\\,c\\right)\\,A_{n}\\,x^{n}}{n!}}\\]"
      ],
      "text/plain": [
       "            3                   2                   2\n",
       "        3  d                2  d                2  d\n",
       "(%o21) x  (--- (F(x))) + c x  (--- (F(x))) + b x  (--- (F(x)))\n",
       "             3                   2                   2\n",
       "           dx                  dx                  dx\n",
       "          2                   2\n",
       "      2  d                2  d                    d\n",
       " + a x  (--- (F(x))) + 3 x  (--- (F(x))) + b c x (-- (F(x)))\n",
       "           2                   2                  dx\n",
       "         dx                  dx\n",
       "          d                 d                   d                 d\n",
       " + a c x (-- (F(x))) + c x (-- (F(x))) + a b x (-- (F(x))) + b x (-- (F(x)))\n",
       "          dx                dx                  dx                dx\n",
       "        d               d\n",
       " + a x (-- (F(x))) + x (-- (F(x))) + a b c F(x) = \n",
       "        dx              dx\n",
       "inf     3                2                                    n\n",
       "====  (n  + (c + b + a) n  + ((b + a) c + a b) n + a b c) A  x\n",
       "\\                                                          n\n",
       " >    ---------------------------------------------------------\n",
       "/                                n!\n",
       "====\n",
       "n = 0"
      ],
      "text/x-maxima": [
       "x^3*'diff(F(x),x,3)+c*x^2*'diff(F(x),x,2)+b*x^2*'diff(F(x),x,2)\n",
       "                   +a*x^2*'diff(F(x),x,2)+3*x^2*'diff(F(x),x,2)\n",
       "                   +b*c*x*'diff(F(x),x,1)+a*c*x*'diff(F(x),x,1)\n",
       "                   +c*x*'diff(F(x),x,1)+a*b*x*'diff(F(x),x,1)\n",
       "                   +b*x*'diff(F(x),x,1)+a*x*'diff(F(x),x,1)+x*'diff(F(x),x,1)\n",
       "                   +a*b*c*F(x)\n",
       "  = 'sum(((n^3+(c+b+a)*n^2+((b+a)*c+a*b)*n+a*b*c)*A[n]*x^n)/n!,n,0,inf)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F14:factor(expand(sumcontract(intosum(FD3_1*x^3+p*FD2_1*x^2+q*FD1_1*x+r*F0)))),SOL2;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "5b815e9a-1062-40e8-9cb9-83f6ff6c709f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\tag{${\\it \\%o}_{22}$}\\left(x^3-x^2\\right)\\,\\left(\\frac{d^3}{d\\,x^3}\\,F\\left(x\\right)\\right)+\\left(\\left(c+b+a+3\\right)\\,x^2+\\left(-e-d-1\\right)\\,x\\right)\\,\\left(\\frac{d^2}{d\\,x^2}\\,F\\left(x\\right)\\right)+\\left(\\left(\\left(b+a+1\\right)\\,c+\\left(a+1\\right)\\,b+a+1\\right)\\,x-d\\,e\\right)\\,\\left(\\frac{d}{d\\,x}\\,F\\left(x\\right)\\right)+a\\,b\\,c\\,F\\left(x\\right)=0\\]"
      ],
      "text/plain": [
       "                   3\n",
       "         3    2   d                               2\n",
       "(%o22) (x  - x ) (--- (F(x))) + ((c + b + a + 3) x  + ((- e) - d - 1) x)\n",
       "                    3\n",
       "                  dx\n",
       "   2\n",
       "  d                                                            d\n",
       " (--- (F(x))) + (((b + a + 1) c + (a + 1) b + a + 1) x - d e) (-- (F(x)))\n",
       "    2                                                          dx\n",
       "  dx\n",
       " + a b c F(x) = 0"
      ],
      "text/x-maxima": [
       "(x^3-x^2)*'diff(F(x),x,3)+((c+b+a+3)*x^2+((-e)-d-1)*x)*'diff(F(x),x,2)\n",
       "                         +(((b+a+1)*c+(a+1)*b+a+1)*x-d*e)*'diff(F(x),x,1)\n",
       "                         +a*b*c*F(x)\n",
       "  = 0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lhs(F14)-lhs(F9)=0,ratsimp;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fca5157b-8930-4b40-90b9-b622a44b851c",
   "metadata": {},
   "source": [
    "求めたかった微分方程式が得られました。\n",
    "\n",
    "より複雑な超幾何関数${}_pF_q$でも$Q(n)\\,A_{n+1}=P(n)\\,A_n$、ただし$P(n)$は$n$の$p$次多項式、$Q(n)$は$n$の$q$次多項式という関係がありますから、$P(n)$を$\\{(n)^-_i | i=0,1,\\dots,p \\}$の要素の1次結合、$Q(n)$を$\\{(n)^-_i | i=0,1,\\dots,q \\}$の要素の1次結合で書くことができて、それらを使って微分の1次結合を求めることができます。\n",
    "\n",
    "なお下降階乗冪が多項式の基底になることはWikipediaの階乗冪の項目にも記載があり参考にしました。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "860fabeb-7ba2-47e1-a55f-8d453e0c1112",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Maxima",
   "language": "maxima",
   "name": "maxima"
  },
  "language_info": {
   "codemirror_mode": "maxima",
   "file_extension": ".mac",
   "mimetype": "text/x-maxima",
   "name": "maxima",
   "pygments_lexer": "maxima",
   "version": "5.44.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}