|
|
!最も簡単な微分方程式 dy/dx=f'(x) の解 f(x)の曲線描画 …不定積分と同、について、
!４次のルンゲクッタ法が、４次のテイラー展開と同精度になる理由。
!------------------------------------------------------------
!区間 h 直接のテイラー展開式
!f(x+h)= f(x) + f'(x)*h + f''(x)/2!*h^2 + f'''(x)/3!*h^3 + f''''(x)/4!*h^4 + ...
!４階微分を要する４次のテイラー展開の出力が、
!区間 h を細分化した位置 (a1,a2,a3,a4) での、１階微分だけの重み付き (b1,b2,b3,b4) 加算でも、
!同じになると仮定し、次の様に置く。
!f(x+h)= f(x)+
! +b1*f'(x+a1*h)*h
! +b2*f'(x+a2*h)*h
! +b3*f'(x+a3*h)*h
! +b4*f'(x+a4*h)*h
!この式を更にテイラー展開し、直接のテイラー展開式と、同じになれるか a1,a2,, b1,b2,, を調べる。
!f(x+h)= f(x)+
! +b1*{f'(x) +f''(x)/1!*(a1*h) +f'''(x)/2!*(a1*h)^2 +f''''(x)/3!*(a1*h)^3 +？*h^4}*h
! +b2*{f'(x) +f''(x)/1!*(a2*h) +f'''(x)/2!*(a2*h)^2 +f''''(x)/3!*(a2*h)^3 +？*h^4}*h
! +b3*{f'(x) +f''(x)/1!*(a3*h) +f'''(x)/2!*(a3*h)^2 +f''''(x)/3!*(a3*h)^3 +？*h^4}*h
! +b4*{f'(x) +f''(x)/1!*(a4*h) +f'''(x)/2!*(a4*h)^2 +f''''(x)/3!*(a4*h)^3 +？*h^4}*h
!各 階微分ごとに係数を整理。
! +(b1+b2+b3+b4) *f'(x)*h
! +(b1*a1 +b2*a2 +b3*a3 +b4*a4 ) *f''(x)*h^2
! +(b1*a1^2+b2*a2^2+b3*a3^2+b4*a4^2)/2*f'''(x)*h^3
! +(b1*a1^3+b2*a2^3+b3*a3^3+b4*a4^3)/6*f''''(x)*h^4
! +？*h^5
!係数比較により、以下の式を満たせば同じになる。
! b1 +b2 +b3 +b4 =1
! (b1*a1 +b2*a2 +b3*a3 +b4*a4 ) =1/2!
! (b1*a1^2+b2*a2^2+b3*a3^2+b4*a4^2)/2=1/3!
! (b1*a1^3+b2*a2^3+b3*a3^3+b4*a4^3)/6=1/4!
! | 1 , 1 , 1 , 1 ||b1| =|1 |
! | a1 , a2 , a3 , a4 ||b2| =|1/2|
! | a1^2, a2^2, a3^2, a4^2||b3| =|1/3|
! | a1^3, a2^3, a3^3, a4^3||b4| =|1/4|
!ルンゲクッタ４は通常、
!a1..a4= 0, 1/2, 1/2, 1 で、逆行列が、求まらないけれども、(特異行列エラー)
!b1..b4= 1/6, 1/3, 1/3, 1/6 は、上式を満たしているのが分る。
! a1,a2,, を適宜に決めると、b1,b2,, が決まる式になっているので、
!逆行列の求まる範囲で、他の a1,a2,, で試しに、グラフの重ね書きをしてみる。
!a1,a2,a3,a4 と、b1,b2,b3,b4 の組合せが無数にあるのが判る。
DIM A_(4,4), I_(4,4), B_(4,1), C_(4,1)
!-----
SET DRAW mode mask
SET COLOR MIX(15) .4,.4,.4
SET WINDOW -2, 9, -5, 3
DRAW grid( EXP(1), 1)
!----
LET xs= EXP(-4) !開始ｘ
LET xe= EXP( 2) !終了ｘ
LET h=.01 !計算間隔
!------------------------------
! 原始関数（比較用） !原始関数、サンプル１、サンプル２
!------------------------------ !３曲線を、減法混色３原色で
SET LINE COLOR "cyan" !重ねるので、一致していれば
FOR x=xs TO xe STEP h !黒線になる。
PLOT LINES: x,f(x);
NEXT x
PLOT LINES
WAIT DELAY .5
!------------------------------
! ルンゲクッタ a1..a4 サンプル１
!------------------------------
CALL a1a4_b1b4( 0, 1/2-.001, 1/2+.001, 1) !a1,a2,a3,a4 から b1,b2,b3,b4 を計算
SET LINE COLOR "yellow"
CALL GRAPH
WAIT DELAY .5
!------------------------------
! ルンゲクッタ a1..a4 サンプル２
!------------------------------
CALL a1a4_b1b4( 0, 1/3, 2/3, 1) !a1,a2,a3,a4 から b1,b2,b3,b4 を計算
SET LINE COLOR "magenta"
CALL GRAPH
SUB a1a4_b1b4( a1,a2,a3,a4 )
CALL A_row(1, 1 , 1 , 1 , 1 )
CALL A_row(2, a1 , a2 , a3 , a4 )
CALL A_row(3, a1^2, a2^2, a3^2, a4^2 )
CALL A_row(4, a1^3, a2^3, a3^3, a4^3 )
MAT I_=INV(A_)
LET C_(1,1)=1
LET C_(2,1)=1/2
LET C_(3,1)=1/3
LET C_(4,1)=1/4
MAT B_=I_*C_
PRINT "b1,b2,b3,b4 の値"
MAT PRINT B_;
END SUB
SUB A_row(r, a,b,c,d)
LET A_(r,1)=a
LET A_(r,2)=b
LET A_(r,3)=c
LET A_(r,4)=d
END SUB
!------------------------------
! 原始関数を使わないで、初期値から、その積分曲線を、ルンゲクッタで、描く。
!------------------------------
SUB GRAPH
LET y=f(xs) !出力初期値
FOR x=xs TO xe STEP h
PLOT LINES: x,y;
CALL RungeKutta
NEXT x
PLOT LINES
END SUB
!-----------------------------
!f(x+h)= f(x)+
! +b1*f'(x+a1*h)*h
! +b2*f'(x+a2*h)*h
! +b3*f'(x+a3*h)*h
! +b4*f'(x+a4*h)*h
!-----------------------------
SUB RungeKutta
CALL Diff( k1, x+A_(2,1)*h)
CALL Diff( k2, x+A_(2,2)*h)
CALL Diff( k3, x+A_(2,3)*h)
CALL Diff( k4, x+A_(2,4)*h)
LET y=y+(B_(1,1)*k1 +B_(2,1)*k2 +B_(3,1)*k3 +B_(4,1)*k4)*h
END SUB
!********* サンプル関数 ********* DEF f()= 、LET k= は、同時交換。
DEF f(x)=LOG(x) !原始関数：f(x)=log(x) １階微分：f'(x)=1/x
!DEF f(x)=SIN(x^2)-COS(x) !原始関数：f(x)=SIN(x^2)-COS(x) １階微分：f'(x)=COS(x^2)*2*x+SIN(x)
!DEF f(x)=x^3-3*x^2+x-1 !原始関数：f(x)=x^3-3*x^2+x-1 １階微分：f'(x)=3*x^2-6*x+1
SUB Diff( k, x)
LET k=1/x !k= dy/dx= f'(x) 最も簡単な微分方程式
!LET k=COS(x^2)*2*x+SIN(x) !k= dy/dx= f'(x) 最も簡単な微分方程式
!LET k=3*x^2-6*x+1 !k= dy/dx= f'(x) 最も簡単な微分方程式
END SUB
END
|
|