Euler And Range-Kutta

• Tweet
• 

from math import floor
from math import cos

def euler(xi, yi, f, xf, h):
    x = xi
    y = yi
    xlist = []
    ylist = []
    
    n = int(floor(xf-xi)/h)

    for i in range(0,n):
        xlist.append(x)
        ylist.append(y)
        m = f(x, y)
        y = y + h*m
        x = x + h
        i+=1
    
    return map(None, xlist, ylist)

def rk(xi, yi, f, xf, h):
        x = xi
        y = yi
        xlist = []
        ylist = []

        n = int(floor(xf-xi)/h)

        for i in range(0,n):
            xlist.append(x)
            ylist.append(y)
            m1 = f(x, y)
            m2 = f(x+h/2, y+m1/2)
            m3 = f(x+h/2, y+m2/2)
            m4 = f(x+h, y+m3)
            y = y + h*(m1/6 + m2/3 + m3/3 + m4/6)
            x = x + h
            i+=1

        return map(None, xlist, ylist)

def f(x, y):
    return (-y+cos(4*x))/4


t = euler(0.0, -1.0, f, 5.0, 0.4)
# # TODO erro
# #
# 
for x,y in t:
     print "%f, %f" % (x,y)


# print "---------------------------"
# # hi = 0.2
# # hf = 0.05
# # h = hi
# # s = []
# # qc = []
# # ec = []
# # 
# # r1 = rk(1.0, 2.0, f, 2.0, 1.0)
# # r2 = rk(1.0, 2.0, f, 2.0, 0.5)
# # r3 = rk(1.0, 2.0, f, 2.0, 0.25)
# # 
# # 
# # qc = (r2[0][1]-r1[0][1]) / (r3[0][1]-r2[0][1])
# # print qc
# 
# # while h >= hf:
# #     r = rk(1.0, 2.0, f, 10.0, h)
# #     s.append(r)
# #     h = h / 2
# # 
# # 
# # for j in range(0, len(s[0])):
# #     for i in range(0, len(s)-2):
# #         div = (s[i+2][4*j][1] - s[i+1][2*j][1]) # S''-S'
# #         if div == 0:
# #             break
# # 
# #         t = (s[i+1][2*j][1] - s[i][j][1]) / div # (S'-S)/div com div = (S''-S')<=>(S'-S)/(S''-S')
# #         qc.append(t)
# #         # t = div / 15
# #         # ec.append(t)
# 
# for i in range(0, len(qc)):
#     print qc[i]
# #    print "%f  ---  %f" % (s[0][i][0], qc[i])
# #qc ~= 16 implica validade => podemos usar os valores do erro: e => qualidade
# # for i in range(0, len(qc)):
# #     print "S=%.12f, S'=%.12f, S''=%.12f\n\t=> qc=%.12f\n" \
# #             % (s[i][], s[i+1], s[i+2], qc[i])

Tags:
• Numerical Methods
• Python