例題 5-11 利用 四階 Runge-Kutta 解 僵硬(Stiff Differential Equation)常微分程式 y' = (50/y) - 50y , 0<= t <=1 , y(0) = 2^0.5 = 1.41414

'''
例題 5-11 利用 四階 Runge-Kutta 解 僵硬(Stiff Differential Equation)常微分程式
  y' =  (50/y) - 50y , 0<= t <=1 , y(0) = 2^0.5 = 1.41414
 
/* Based on Four-Order Runge-Kutta
 * Method to approximate the solution of the
 * initial-value problem
 *   y'=f(y,t), a<=t<=b, y(a)=y0
 * at (n+1) equally spaced numbers in the interval
 * [a,b]: input a,b,n,and initial condition y0.
 */
'''
import math
def F(y,t):
    return ((50.0/y)-50*y)

def W(t):
    return (math.sqrt(1+math.exp(-100*t) ) )

#======== main========
n1=[1000,100,20,10]
for j in range (0,4):
    n=n1[j]

    a=0.0
    b=1.0
    y0=math.sqrt(2)
    h=(b-a)/n
    t=a
    y=y0
    err=abs(y-W(t))
    print("t \t\ty(t) \t\t    w(t) \t\t   error");
    print("====================================================");
    print("{%.2f} \t\t{%10.7f} \t\t{%10.7f} \t\t{%10.7f} " %(t,y,W(t),err) )
    for i in range (1,n+1):
        k1=h*F(y,t)
        k2=h*F((y+k1/2.0),(t+h/2.0))
        k3=h*F((y+k2/2.0),(t+h/2.0))
        k4=h*F((y+k3),(t+h))
        y=y+(k1+2*k2+2*k3+k4)/6.0
        t=a+i*h
        err=abs(y-W(t))
        if(i% (n/10)==0):
            print("{%.2f}\t\t {%10.7f} \t\t {%10.7f} \t\t {%10.7f} " %(t,y,W(t),err) )

====== RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX5-11.py ========
t y(t)     w(t)    error
====================================================
{0.00} { 1.4142136} { 1.4142136} { 0.0000000} 
{0.10} { 1.0000227} { 1.0000227} { 0.0000000} 
{0.20} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.30} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.40} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.50} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.60} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.70} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.80} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.90} { 1.0000000} { 1.0000000} { 0.0000000} 
{1.00} { 1.0000000} { 1.0000000} { 0.0000000} 
t y(t)     w(t)    error
====================================================
{0.00} { 1.4142136} { 1.4142136} { 0.0000000} 
{0.10} { 1.0000270} { 1.0000227} { 0.0000043} 
{0.20} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.30} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.40} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.50} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.60} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.70} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.80} { 1.0000000} { 1.0000000} { 0.0000000} 
{0.90} { 1.0000000} { 1.0000000} { 0.0000000} 
{1.00} { 1.0000000} { 1.0000000} { 0.0000000} 
t y(t)     w(t)    error
====================================================
{0.00} { 1.4142136} { 1.4142136} { 0.0000000} 
{0.10} {31.3406218} { 1.0000227} {30.3405991} 
{0.20} {12.7757104} { 1.0000000} {11.7757104} 
{0.30} { 4.2971794} { 1.0000000} { 3.2971794} 
{0.40} { 1.2416323} { 1.0000000} { 0.2416323} 
{0.50} { 1.9840741} { 1.0000000} { 0.9840741} 
{0.60} { 5.1231349} { 1.0000000} { 4.1231349} 
{0.70} { 8.9870724} { 1.0000000} { 7.9870724} 
{0.80} { 2.0245028} { 1.0000000} { 1.0245028} 
{0.90} { 6.3456384} { 1.0000000} { 5.3456384} 
{1.00} {-1.9624394} { 1.0000000} { 2.9624394} 
t y(t)     w(t)    error
====================================================
{0.00} { 1.4142136} { 1.4142136} { 0.0000000} 
{0.10} {-15.8525966} { 1.0000227} {16.8526193} 
{0.20} {-215.7458639} { 1.0000000} {216.7458639} 
{0.30} {-2957.4012711} { 1.0000000} {2958.4012711} 
{0.40} {-40541.0340389} { 1.0000000} {40542.0340389} 
{0.50} {-555750.0076712} { 1.0000000} {555751.0076712} 
{0.60} {-7618406.3551152} { 1.0000000} {7618407.3551152} 
{0.70} {-104435653.7847015} { 1.0000000} {104435654.7847015} 
{0.80} {-1431638753.9652824} { 1.0000000} {1431638754.9652824} 
{0.90} {-19625381252.2740898} { 1.0000000} {19625381253.2740898} 
{1.00} {-269031267999.9240723} { 1.0000000} {269031268000.9240723} 
>>>