2019年1月23日 星期三

例題 5-6 利用 四階 Runge-Kutta 解一階常微分方程式 y' = - y + t^2 +1 , y(0)=1 , 0<= t <= 1 取 h=0.01 真實解 W(t) = -2 exp(-t) + t^2 - 2t +3

'''
#例題 5-6 利用 四階 Runge-Kutta  解一階常微分方程式
# y' = - y + t^2  +1  , y(0)=1  , 0<= t <= 1
# 取 h=0.01
# 真實解 W(t) =  -2 exp(-t)  + t^2 - 2t +3

誤差 err = { 0.0000000} 

    Based on Four-Order Runge-Kutta  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+t*t+1)
def W(t):
    return (-2*(1.0/math.exp(t))+pow(t,2)-2*t+3)

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




======== RESTART: F:/2018-09勤益科大數值分析/數值分析/PYTHON/EX5-6.py ============
      t               y(t)                    w(t)             error
=========================================================
{0.00}   { 1.0000000}    { 1.0000000} { 0.0000000} 
{0.10}   { 1.0003252}    { 1.0003252} { 0.0000000} 
{0.20}   { 1.0025385}    { 1.0025385} { 0.0000000} 
{0.30}   { 1.0083636}    { 1.0083636} { 0.0000000} 
{0.40}   { 1.0193599}    { 1.0193599} { 0.0000000} 
{0.50}   { 1.0369387}    { 1.0369387} { 0.0000000} 
{0.60}   { 1.0623767}    { 1.0623767} { 0.0000000} 
{0.70}   { 1.0968294}    { 1.0968294} { 0.0000000} 
{0.80}   { 1.1413421}    { 1.1413421} { 0.0000000} 
{0.90}   { 1.1968607}    { 1.1968607} { 0.0000000} 
{1.00}   { 1.2642411}    { 1.2642411} { 0.0000000} 
>>> 

沒有留言:

張貼留言

Messaging API作為替代方案

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