例題 5-3 利用修正 Euler 法 解 一階常微分方程式 y' - y = 2 exp(4t) , y(0)=-3 , 0<= t <= 1 取 h=0.01 真實解 W(t) = (2/3) exp(4t) - (11/3) exp(t)

Modified Euler's Method :

The Euler forward scheme may be very easy to implement but it can't give accurate solutions.    A  very small step size is required for any meaningful result.  In this scheme, since, the starting point of each sub-interval is used to find the slope of the solution curve,  the solution would be correct only if the function is linear. So an improvement over this is to take the arithmetic average of the slopes at x_i  and x_i+1(that is, at the end points of each sub-interval). The scheme so obtained is called modified Euler's method. It works first by approximating a value to y_i+1 and then improving it by making use of average slope.

yi+1	= yi+ h/2 (y'i + y'i+1)
	= yi + h/2(f(xi, yi) + f(xi+1, yi+1))

 If Euler's method is used to find the first approximation of y_i+1 then 

y_i+1 = y_i + 0.5h(f_i  + f(x_i+1, y_i + hf_i))
源自於
https://mat.iitm.ac.in/home/sryedida/public_html/caimna/ode/euler/ie.html

'''
#例題 5-3 利用修正 Euler 法 解 一階常微分方程式
# y' - y = 2 exp(4t) , y(0)=-3  , 0<= t <= 1
# 取 h=0.01
# 真實解 W(t) = (2/3) exp(4t) - (11/3) exp(t)
/* Based on modified Euler
 * 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 (y+ 2*math.exp(4*t) )

def w(t):
    return ((2.0/3.0)*math.exp(4.0*t)-(11.0/3)*math.exp(t))

#======== main========
n=100
a=0.0
b=1.0
y0=-3.0
h=(b-a)/n
t=a
y=y0
print("t \t\t y(t) \t\t\t       w(t) \t\t\t  error");
print("=========================================================")
print("{%.2f} \t\t  {%10.6f}  \t\t  {%10.6f} \t\t {%10.6f} " %(t,y,w(t),abs(y-w(t)) ) )
for i in range (1,n+1):
    k1=h*F(y,t)
    k2=h*F((y+k1),(t+h))
    y=y+(k1+k2)/2.0
    t=a+i*h;
    if(i%10==0):
       print("{%.2f} \t\t  {%10.6f}  \t\t  {%10.6f} \t\t {%10.6f} " %(t,y,w(t),abs(y-w(t)) ) )

======== RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX5-3.py =============
t y(t)        w(t)   error
=========================================================
{0.00}   { -3.000000}    { -3.000000} {  0.000000} 
{0.10}   { -3.057725}    { -3.057744} {  0.000018} 
{0.20}   { -2.994738}    { -2.994783} {  0.000045} 
{0.30}   { -2.735987}    { -2.736071} {  0.000084} 
{0.40}   { -2.167862}    { -2.168002} {  0.000140} 
{0.50}   { -1.119051}    { -1.119274} {  0.000223} 
{0.60}   {  0.668025}    {  0.667682} {  0.000343} 
{0.70}   {  3.579857}    {  3.579338} {  0.000519} 
{0.80}   {  8.195482}    {  8.194703} {  0.000779} 
{0.90}   { 15.381439}    { 15.380278} {  0.001161} 
{1.00}   { 26.433458}    { 26.431733} {  0.001725} 
>>>