2019年1月23日 星期三

例題 5-2 利用向前 Euler 解 一階常微分方程式 y' - y = 2 exp(4t) , y(0)=-3 , 0<= t <= 1 取 h=0.01 , 0.001 , 0.0001

'''
#例題 5-2 利用向前 Euler 解 一階常微分方程式
# y' - y = 2 exp(4t) , y(0)=-3  , 0<= t <= 1
# 取 h=0.01 , 0.001 , 0.0001
# 真實解 W(t) = (2/3) exp(4t) - (11/3) exp(t)
 * Forward Euler Method is used for
 * solving y'=f(y,t) of first order  ordinary differential equation
 * with initial condition y(0)=y0  known.
 *
 '''
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========
t0=0
t1=1
h=[0.01, 0.001 ,0.0001]
n=[0,0,0]
m=[0,0,0]
y0=-3
for i in range (0,3):
    n[i]=int ( (t1-t0)/h[i] )
    m[i]=int(n[i]/10)
#print (n,m)

for j in range (0,3):   
    y=y0;
    t=t0;

    print("t \t\t y(t) \t\t\t       w(t) \t\t\t  error");
    print("=========================================================")
    print('h=',h[j])
    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[j]+1):
        t=t+h[j]
        y=y+h[j]*F(y,t)
        if (i %m[j]==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-2.py ============
t y(t)        w(t)   error
=========================================================
h= 0.01
{0.00}   { -3.000000}    { -3.000000} {  0.000000} 
{0.10}   { -3.052267}    { -3.057744} {  0.005477} 
{0.20}   { -2.981340}    { -2.994783} {  0.013443} 
{0.30}   { -2.711053}    { -2.736071} {  0.025018} 
{0.40}   { -2.126148}    { -2.168002} {  0.041855} 
{0.50}   { -1.052877}    { -1.119274} {  0.066397} 
{0.60}   {  0.769944}    {  0.667682} {  0.102262} 
{0.70}   {  3.734158}    {  3.579338} {  0.154820} 
{0.80}   {  8.426750}    {  8.194703} {  0.232046} 
{0.90}   { 15.726079}    { 15.380278} {  0.345801} 
{1.00}   { 26.945465}    { 26.431733} {  0.513732} 
t y(t)        w(t)   error
=========================================================
h= 0.001
{0.00}   { -3.000000}    { -3.000000} {  0.000000} 
{0.10}   { -3.057197}    { -3.057744} {  0.000546} 
{0.20}   { -2.993442}    { -2.994783} {  0.001341} 
{0.30}   { -2.733576}    { -2.736071} {  0.002495} 
{0.40}   { -2.163830}    { -2.168002} {  0.004172} 
{0.50}   { -1.112658}    { -1.119274} {  0.006616} 
{0.60}   {  0.677869}    {  0.667682} {  0.010187} 
{0.70}   {  3.594757}    {  3.579338} {  0.015419} 
{0.80}   {  8.217807}    {  8.194703} {  0.023104} 
{0.90}   { 15.414700}    { 15.380278} {  0.034421} 
{1.00}   { 26.482860}    { 26.431733} {  0.051126} 
t y(t)        w(t)   error
=========================================================
h= 0.0001
{0.00}   { -3.000000}    { -3.000000} {  0.000000} 
{0.10}   { -3.057689}    { -3.057744} {  0.000055} 
{0.20}   { -2.994649}    { -2.994783} {  0.000134} 
{0.30}   { -2.735822}    { -2.736071} {  0.000249} 
{0.40}   { -2.167585}    { -2.168002} {  0.000417} 
{0.50}   { -1.118613}    { -1.119274} {  0.000661} 
{0.60}   {  0.668700}    {  0.667682} {  0.001018} 
{0.70}   {  3.580879}    {  3.579338} {  0.001541} 
{0.80}   {  8.197013}    {  8.194703} {  0.002309} 
{0.90}   { 15.383719}    { 15.380278} {  0.003441} 
{1.00}   { 26.436843}    { 26.431733} {  0.005110} 
>>> 

沒有留言:

張貼留言

Messaging API作為替代方案

  LINE超好用功能要沒了!LINE Notify明年3月底終止服務,有什麼替代方案? LINE Notify將於2025年3月31日結束服務,官方建議改用Messaging API作為替代方案。 //CHANNEL_ACCESS_TOKEN = 'Messaging ...