s2Eva_IT2012_T2 Modelo de clima

Se deja de tarea realizar las tres primeras iteraciones en papel.

Se presenta el resultado usando el algoritmo de segundo grado como una variante a la respuesta usada como ejemplo.

 [ ti, xi, yi, zi]
[[ 0.0000e+00  1.0000e+01  7.0000e+00  7.0000e+00]
 [ 2.5000e-03  9.9323e+00  7.5033e+00  7.1335e+00]
 [ 5.0000e-03  9.8786e+00  7.9988e+00  7.2774e+00]
 ...
 [ 2.4995e+01 -8.4276e+00 -2.7491e+00  3.3021e+01]
 [ 2.4998e+01 -8.2860e+00 -2.6392e+00  3.2858e+01]
 [ 2.5000e+01 -8.1453e+00 -2.5346e+00  3.2692e+01]]

Algoritmo en Python

 2Eva_IT2012_T2 Modelo de clima
# MATG1013 Análisis Numérico
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,j,t0,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [t0,x0,y0,z0]
    estimado[0] = [t0,x0,y0,z0]
    ti = t0
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi,zi)
        K1y = h * g(ti,xi,yi,zi)
        K1z = h * j(ti,xi,yi,zi)
        
        K2x = h * f(ti+h,xi + K1x, yi + K1y, zi + K1z)
        K2y = h * g(ti+h,xi + K1x, yi + K1y, zi + K1z)
        K2z = h * j(ti+h,xi + K1x, yi + K1y, zi + K1z)
        
        xi = xi + (K1x+K2x)/2
        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        ti = ti + h
        
        estimado[i] = [ti,xi,yi,zi]
    return(estimado)


#INGRESO
to = 0
xo = 10
yo = 7
zo = 7
f = lambda t,x,y,z: 10*(y-x)
g = lambda t,x,y,z: x*(28-z) - y
j = lambda t,x,y,z: -(8/3)*z + (x*y)
h = 0.0025
muestras = 10000

#PROCEDIMIENTO
#Rugen-Kutta 2_orden
tabla = rungekutta2_fg(f,g,j,to,xo,yo,zo,h,muestras)

#SALIDA
np.set_printoptions(precision=4)
print(' [ ti, xi, yi, zi]')
print(tabla)

# Gráfica
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(tabla[:,1], tabla[:,2], tabla[:,3])
plt.show()