例題 5-10 利用 四階 Runge-Kutta 解 三個聯立常微分程式 y1' = y2 - y3 + t , y1(0)=1 y2' = 3t^2 , y2(0)=1 , 0<= t <=1 y3' = y2 + esp (-t) , y3(0)=-1

'''
例題 5-10 利用 四階 Runge-Kutta 解 三個聯立常微分程式
     y1' = y2 - y3 + t     ,  y1(0)=1
     y2' = 3t^2               ,  y2(0)=1   ,  0<=  t  <=1
     y3' = y2 + esp (-t)  ,  y3(0)=-1
真實解
W1(t) = -0.05 t^5 + 0.25t^4 + t + 2 - exp(-t)
W2(t) = t^3 +1
W3(t) = 0.25 t^4 + t - exp(-t)

0<=  t  <=1 計算 y1(x) , y2(x) , y3(t) 設 h=0.1 , 0.01
   
/* ex5-10.c based on Fourth-Order Runge-Kutta Method
 * to approximate the solution of the m=3 order
 * system of first-order initial-value problem.
 *    y1=f1(y1,y2,...,ym,t)
 *    y2=f2(y1,y2,...,ym,t)
 *        .
 *        .
 *    ym=fm(y1,y2,...,ym,t)
 *    a<=t<=b, y1(a)=y01,y2(a)=y02,...,ym(a)=y0m.
 * at (n+1) equally spaced numbers in the interval
 * [a,b].
 */
 '''
import math
def f1(y1,y2,y3,t):
    return (y2-y3+t)

def f2(y1,y2,y3,t):
    return (3*pow(t,2))

def f3(y1,y2,y3,t):
    return (y2+(1.0/math.exp(t)))

def w1(t):
    return (-0.05*pow(t,5)+0.25*pow(t,4)+t+2.0-(1.0/math.exp(t)))

def w2(t):
    return (pow(t,3)+1.0)

def w3(t):
    return (0.25*pow(t,4)+t-(1.0/math.exp(t)))

#======== main========
n=10
m=3
a=0.0
b=1.0
k1=[  0.0 for i in range (11) ]
k2=[  0.0 for i in range (11) ]
k3=[  0.0 for i in range (11) ]
k4=[  0.0 for i in range (11) ]
y=[  0.0 for i in range (11) ]
y0=[  0.0 for i in range (11) ]
h=(b-a)/n
t=a
y0[1]=1.0
y0[2]=1.0
y0[3]=-1.0
for j in range (1 , m +1 ):
    y[j]=y0[j]

print("t \t   y1 \t   |y1-w1(t)|  \ty2  \t  |y2-w2(t)| \t y3  \t  |y3-w3(t)|");
for j in range (0,78):
    print('=',end='')
print() 
print("{%.2f} \t  {%10.7f}  \t {%10.7f} \t {%10.7f}  \t {%10.7f}\t {%10.7f}\t {%10.7f}" %(t,y0[1],abs(y0[1]-w1(t)),y0[2],abs(y0[2]-w2(t)) ,y0[3],abs(y0[3]-w3(t)) ))
for i in range (1,n+1):
    for j in range(1 , m+1):
        if(j==1):
            k1[j]=h*f1((y[j]),(y[j+1]),(y[j+2]),t);
        elif(j==2):
            k1[j]=h*f2((y[j-1]),(y[j]),(y[j+1]),t);
        elif(j==3):
            k1[j]=h*f3((y[j-2]),(y[j-1]),(y[j]),t);

    for j in range (1 , m+1):
        if(j==1):
            k2[j]=h*f1((y[j]+0.5*k1[j]),(y[j+1]+0.5*k1[j+1]),(y[j+2]+0.5*k1[j+2]),(t+0.5*h))
        elif(j==2):
            k2[j]=h*f2((y[j-1]+0.5*k1[j-1]),(y[j]+0.5*k1[j]),(y[j+1]+0.5*k1[j+1]),(t+0.5*h))
        elif(j==3):
            k2[j]=h*f3((y[j-2]+0.5*k1[j-2]),(y[j-1]+0.5*k1[j-1]),(y[j]+0.5*k1[j]),(t+0.5*h));

    for j in range (1 , m+1):
        if(j==1):
            k3[j]=h*f1((y[j]+0.5*k2[j]),(y[j+1]+0.5*k2[j+1]),(y[j+2]+0.5*k2[j+2]),(t+0.5*h))
        elif(j==2):
            k3[j]=h*f2((y[j-1]+0.5*k2[j-1]),(y[j]+0.5*k2[j]),(y[j+1]+0.5*k2[j+1]),(t+0.5*h))
        elif(j==3):
            k3[j]=h*f3((y[j-2]+0.5*k2[j-2]),(y[j-1]+0.5*k2[j-1]),(y[j]+0.5*k2[j]),(t+0.5*h))

    for j in range (1 , m+1):
        if(j==1):
            k4[j]=h*f1((y[j]+k3[j]),(y[j+1]+k3[j+1]),(y[j+2]+k3[j+2]),(t+h))
        elif(j==2):
            k4[j]=h*f2((y[j-1]+k3[j-1]),(y[j]+k3[j]),(y[j+1]+k3[j+1]),(t+h))
        elif(j==3):
            k4[j]=h*f3((y[j-2]+k3[j-2]),(y[j-1]+k3[j-1]),(y[j]+k3[j]),(t+h))

    for j in range (1 , m+1):
        y[j]=y[j]+((k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/6.0)
    t=a+i*h;
    if(i%1==0):
        print("{%.2f} \t  {%10.7f}  \t {%10.7f} \t {%10.7f}  \t {%10.7f}\t {%10.7f}\t {%10.7f}" %(t,y[1],abs(y[1]-w1(t)),y[2],abs(y[2]-w2(t)) ,y[3],abs(y[3]-w3(t)) ))




#======== main========
n=10

    if(i%1==0):

#======== main========
n=100

    if(i%10==0):

========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX5-10.py ===========
t    y1    |y1-w1(t)|  y2    |y2-w2(t)| y3    |y3-w3(t)|
=========================================================
{0.00} { 1.0000000} { 0.0000000} { 1.0000000} { 0.0000000} {-1.0000000} { 0.0000000}
{0.10} { 1.1951869} { 0.0000001} { 1.0010000} { 0.0000000} {-0.8048124} { 0.0000000}
{0.20} { 1.3816530} { 0.0000003} { 1.0080000} { 0.0000000}{-0.6183307} { 0.0000000}
{0.30} { 1.5610849} { 0.0000004} { 1.0270000} { 0.0000000}{-0.4387932} { 0.0000000}
{0.40} { 1.7355674} { 0.0000005} { 1.0640000} { 0.0000000}{-0.2639200} { 0.0000000}
{0.50} { 1.9075312} { 0.0000007} { 1.1250000}  { 0.0000000}{-0.0909056} { 0.0000000}
{0.60} { 2.0796995} { 0.0000008} { 1.2160000} { 0.0000000}{ 0.0835884} { 0.0000000}
{0.70} { 2.2550352} { 0.0000010} { 1.3430000}  { 0.0000000}{ 0.2634397} { 0.0000000}
{0.80} { 2.4366860} { 0.0000011} { 1.5120000}  { 0.0000000}{ 0.4530711} { 0.0000000}
{0.90} { 2.6279296} { 0.0000012} { 1.7290000} { 0.0000000}{ 0.6574554} { 0.0000000}
{1.00} { 2.8321192} { 0.0000014} { 2.0000000}  { 0.0000000}{ 0.8821206} { 0.0000000}
>>> 

#======== main========
n=100

    if(i%10==0):

========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX5-10.py ===========
t    y1    |y1-w1(t)|  y2    |y2-w2(t)| y3    |y3-w3(t)|
===========================================================================
{0.00}   { 1.0000000}  { 0.0000000} { 1.0000000}  { 0.0000000} {-1.0000000} { 0.0000000}
{0.10}   { 1.1951871}  { 0.0000000} { 1.0010000}  { 0.0000000} {-0.8048124} { 0.0000000}
{0.20}   { 1.3816532}  { 0.0000000} { 1.0080000}  { 0.0000000} {-0.6183308} { 0.0000000}
{0.30}   { 1.5610853}  { 0.0000000} { 1.0270000}  { 0.0000000} {-0.4387932} { 0.0000000}
{0.40}   { 1.7355680}  { 0.0000000} { 1.0640000}  { 0.0000000} {-0.2639200} { 0.0000000}
{0.50}   { 1.9075318}  { 0.0000000} { 1.1250000}  { 0.0000000} {-0.0909057} { 0.0000000}
{0.60}   { 2.0797004}  { 0.0000000} { 1.2160000}  { 0.0000000} { 0.0835884} { 0.0000000}
{0.70}   { 2.2550362}  { 0.0000000} { 1.3430000}  { 0.0000000} { 0.2634397} { 0.0000000}
{0.80}   { 2.4366870}  { 0.0000000} { 1.5120000}  { 0.0000000} { 0.4530710} { 0.0000000}
{0.90}   { 2.6279308}  { 0.0000000} { 1.7290000}  { 0.0000000} { 0.6574553} { 0.0000000}
{1.00}   { 2.8321206}  { 0.0000000} { 2.0000000}  { 0.0000000} { 0.8821206} { 0.0000000}

>>>