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

'''
例題 5-7 , 5-8 利用 二階 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-8.c based on Second-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) ]
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.;
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(0,y[j+1],0)
        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(0,(y[j+1]+k1[j+1]),0)
        elif(j==2):
            k2[j]=h*f2((y[j-1]+k1[j-1]),(y[j]+k1[j]),(t+h))
     
    for j in range (1 , m+1):
        y[j]=y[j]+0.5*(k1[j]+k2[j]);
    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-7-5-8.py ===========
t               y1    |y1-w1(t)|         y2            |y2-w2(t)|
============================================================================
{0.00}   { 1.0000000}    { 0.0000000}  { 0.0000000}  { 0.0000000}
{0.10}   { 1.0000000}    { 0.0004372}  { 0.0139758}   { 0.0014486}
{0.20}   { 1.0025157}    { 0.0005377}  { 0.0434233}  { 0.0018048}
{0.30}   { 1.0085206}    { 0.0004548} { 0.0787379}  { 0.0015000}
{0.40}   { 1.0181887}    { 0.0002902} { 0.1131026}  { 0.0008337}
{0.50}   { 1.0311357}    { 0.0001090}  { 0.1418326}  { 0.0000089}
{0.60}   { 1.0466026}    { 0.0000494}  { 0.1618678}  { 0.0008405}
{0.70}   { 1.0635963}    { 0.0001638} { 0.1713791}  { 0.0016304}
{0.80}   { 1.0809948}    { 0.0002250}  { 0.1694648}  { 0.0023110}
{0.90}   { 1.0976266}    { 0.0002318}  { 0.1559163}  { 0.0028568}
{1.00}   { 1.1123295}    { 0.0001879} { 0.1310380}  { 0.0032585}
>>>