sn、cn、dn関数の値を計算する [数値解析]
sn、cn、dn関数の値を計算する
とし、この逆関数
とし、
とおくと、この導関数は
となる。
ここで、
であり、
である。
これらの関係から、sn(u)について次の微分方程式を得ることができる。
ここで改めて、x=u、y=sn(u)とおくと、
という微分方程式を得る。
x=0のとき、
であるので、オイラー法やルンゲ・クッタ法などを用いて数値積分することにより、y=sn(x)の数値解を求めることができる。
具体的には、
とおき、次の連立微分方程式
を、初期条件
のもとで数値的に解けばよい。
k=0.8について、解いた結果は、次の通り。
計算に使用したFortranプログラムは次のとおり。
parameter(n=10)
real x(0:n) , y(0:n), p(0:n)
real k
k=0.8
h=0.1
! 配列変数の初期化
x=0; y=0; p=0
! 計算点の座標を設定
do i= 0, n
x(i)=i*h
end do
! 初期値を設定
y(0) = 0.; p(0)=1.0
! ルンゲ・クッタ法で連立微分方程式を解く
call RK(x,y,p,n)
! 計算結果の出力
write(*,*) ' x sn(x) cn(x) dn(x)'
do i= 0 , n ,1
dn=sqrt(1.0-(k*y(i))**2)
cn=p(i)/dn
write(*,100) x(i),y(i) , cn , dn
end do
100 format(f8.5, 3(1x , f8.5))
end
! ルンゲ・クッタ法(4次)
subroutine RK(x , y, z ,n)
real x(0:n) , y(0:n), z(0:n)
do i= 1 , n
h = x(i)-x(i-1)
dy1 = f(x(i-1),y(i-1),z(i-1))*h
dz1= g(x(i-1),y(i-1),z(i-1))*h
dy2 = h*f(x(i-1)+h/2, y(i-1)+dy1/2, z(i-1)+dz1/2)
dz2 = h*g(x(i-1)+h/2, y(i-1)+dy1/2, z(i-1)+dz1/2)
dy3 = h*f(x(i-1)+h/2, y(i-1)+dy2/2, z(i-1)+dz2/2)
dz3 = h*g(x(i-1)+h/2, y(i-1)+dy2/2, z(i-1)+dz2/2)
dy4 = h*f(x(i-1)+h, y(i-1)+dy3, z(i-1)+dz3)
dz4 = h*g(x(i-1)+h, y(i-1)+dy3, z(i-1)+dz3)
y(i)=y(i-1)+ (dy1 + 2*dy2 + 2*dy3 + dy4)/6.
z(i)=z(i-1)+ (dz1 + 2*dz2 + 2*dz3 + dz4)/6.
end do
end
! dy/dx=f(x,y,z)のf(x,y,z)を定義
function f(x,y,z)
f=z
end
! dz/dx = g(x,y,z)のg(x,y,z)を定義
function g(x, y, z)
real k
k=0.8
g=-(1+k**2)*y + 2.*k**2*y**3
end
4次のルンゲ・クッタ法を使ってsn(x)を数値的に解いたあと、
を計算したあと、
を用いてcs関数を求めてあるにゃ。
少し無駄が多いので、もう少し綺麗にプログラムを書くこともできるだろうけれど(^^ゞ
k=0.8のときの計算結果
CASIOの高精度計算サイトの計算結果
なお、この計算は次のようにしてCASIOさんのところで計算したもの。
CASIOさんがどのようにしてこれを計算しているのかは知らないけれど、このように他人(ひと)の力を頼らずに計算できるというわけ。
コメント 0