例題 5-9 利用 四階 Runge-Kutta 解 二階常微分程式 y'' + 4y' + 4y = 4 cos(t) + 3 sin(t) y(0)=1 , y'(0)=0

'''
例題 5-9 利用 四階 Runge-Kutta 解 二階常微分程式
y'' + 4y' + 4y = 4 cos(t) + 3 sin(t)
y(0)=1  , y'(0)=0

設 z= y'
z' + 4z + 4y = 4 cos(t) + 3 sin(t)
z ' = f2( y , z , t) = - 4z - 4y = 4 cos(t) + 3 sin(t)
f' = f1(y,z ,t) = z  ,  y(0)=1

設 y1= y , y2=z
y'1 = f1(y1,y2,t) = y2  , y1(0)=1
y'2 = f2(y1,y2,t) = - 4y1 - 4y2 + 4 cos(t) + 3 sin(t)  , y2(0)=0

/* ex5-9.c based on Fourth-Order Runge-Kutta Method
 * to approximate the solution of the m 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,t):
    return (y2)

def f2(y1,y2,t):
    return ((-4.0*y2-4.0*y1)+4*math.cos(t)+3*math.sin(t))

def w1(t):
    return ((1+t)*(1.0/math.exp(2*t))+math.sin(t))

def w2(t):
    return (math.cos(t)-(1.0/math.exp(2*t))*(2*t+1))

#======== main========
n=10
m=2;
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]=0.0
for j in range (1 , m +1 ):
    y[j]=y0[j]

print("t \t\t      y1 \t\t\t   |y1-w1(t)|  \t\t\ty2  \t\t  |y2-w2(t)|");
for j in range (0,76):
    print('=',end='')
print() 
print("{%.2f} \t\t  {%10.7f}  \t\t  {%10.7f} \t\t {%10.7f}  \t\t {%10.7f}" %(t,y0[1],abs(y0[1]-w1(t)),y0[2],abs(y0[2]-w2(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],t)
        elif(j==2):
            k1[j]=h*f2(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]),(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]),(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]),(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]),(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]),(t+h))
        elif(j==2):
            k4[j]=h*f2((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\t  {%10.7f}  \t\t  {%10.7f} \t\t {%10.7f}  \t\t {%10.7f}" %(t,y[1],abs(y[1]-w1(t)),y[2],abs(y[2]-w2(t))))

========= RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX5-9.py ============
t                   y1               |y1-w1(t)|  y2    |y2-w2(t)|
=========================================================
{0.00}   { 1.0000000}    { 0.0000000}  { 0.0000000}   { 0.0000000}
{0.10}   { 1.0004348}    { 0.0000024} { 0.0125336}   { 0.0000063}
{0.20}   { 1.0030501}    { 0.0000032} { 0.0416271}   { 0.0000086}
{0.30}   { 1.0089723}    { 0.0000030}  { 0.0772460}   { 0.0000082}
{0.40}   { 1.0184767}    { 0.0000022} { 0.1122749}   { 0.0000061}
{0.50}   { 1.0312436}    { 0.0000011}  { 0.1418267}   { 0.0000031}
{0.60}   { 1.0465534}    { 0.0000002}  { 0.1627080}   { 0.0000004}
{0.70}   { 1.0634340}    { 0.0000015} { 0.1730056}   { 0.0000039}
{0.80}   { 1.0807725}    { 0.0000027} { 0.1717685}   { 0.0000073}
{0.90}   { 1.0973986}    { 0.0000038}  { 0.1587628}   { 0.0000103}
{1.00}   { 1.1121462}    { 0.0000047} { 0.1342835} { 0.0000130}
>>>