ロピタルの定理 [fortran]

こんにちは、管理人のさくらえびです。

ロピタルの定理をプログラムでやりたいと思います。
ロピタルの定理といえば、大学受験での使用について議論されていますが、実は、17-18世紀ごろから、ベルヌーイとロピタルのどちらが作ったかという話で論争しております。
何かと論争の多い定理ですね。
しかし、論争が多いということは、非常に有用な定理であるということです。
0近傍での収束値を求めるという点では、
数値計算では重宝します。特に0近傍を扱う計算では、しばしば用いられます。
まずは定義式。cを含む開区間における微分可能な関数\(f(x)\),\(g(x)\)に対し、
\(\lim_{x \to 0} f(c)\),\(\lim_{x \to 0} g(c)\)が互いに０または、±∞であり、かつ
\[
\lim_{x \to 0} f'(x)/g'(x)
\]
が存在し、かつ\(g’(x)’\)≠0ならば、以下の等式が成り立つことを主張する。
\[
\lim_{x \to 0} f(x)/g(x) = \lim_{x \to 0} f'(x)/g'(x)
\]
0点近傍処理が楽になる定理なので、非常に重宝します。
では、さっそくプログラム

function L_HOP(func1,func2)
implicit none
external func1,func2
double precision L_HOP,func1,func2
double precision DIF_F
double precision df1,df2
df1 = DIF_F(func1,0.d0)
df2 = DIF_F(func2,0.d0)
L_HOP = df1/df2
end function

実際に、\(f(x)\)と\(g(x)\)を以下のようにおいてみる。
\[
f(x)=sin(2x)
\]
\[
g(x)=x + sin(x)
\]
このとき、期待する値はロピタルの定理によれば、1なので、
下記のようにテストコードを組んでみます。

program test
implicit none
external fx,gx
double precision x,DIF_F,L_HOP
double precision fx,gx
x = 0.d0
print *, "x,f'(0)/g'(0),f(0)/g(0) =",x,L_HOP(fx,gx),fx(0.d0)/gx(0.d0)
end program
function fx(x)
implicit none
double precision fx,x
fx = sin(2*x)
end function
function gx(x)
implicit none
double precision gx,x
gx = x + sin(x)
end function

結果はこちら

\(f'(x)/g'(x)\)は1に収束し、\(f(x)/g(x)\)は非数になりました。いわゆる、0/0となります。
また、便利な関数を作っていきます。また、C言語でも同様に作っていくので、見ていってください。

by さくらえび