6.2.2 EDO Runge-Kutta d2y/dx2 con Python

[ Runge Kutta d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]
..

1. EDO Runge-Kutta para Segunda derivada d²y/dx²

Para una ecuación diferencial de segunda derivada (segundo orden) con condiciones de inicio en x₀, y₀, y’₀

\frac{\delta ^2 y}{\delta x^2} = \frac{\delta y}{\delta x} + etc

Forma estandarizada de la ecuación:

y'' = y' + etc

Se puede sustituir la variable y’ por z, lo que se convierte a dos expresiones que forman un sistema de ecuaciones:

\begin{cases} z= y' = f_x(x,y,z) \\ z' = (y')' = z + etc = g_x(x,y,z) \end{cases}

y se pueden reutilizar los métodos para primeras derivadas, por ejemplo Runge-Kutta de 2do y 4to orden para las variables x,y,z de forma simultanea.

Runge-Kutta 2do Orden tiene error de truncamiento O(h³)
Runge-Kutta 4do Orden tiene error de truncamiento O(h⁵)

[ Runge Kutta d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]

2. Runge-Kutta 2do Orden para Segunda derivada d²y/dx² en Python

y'' = y' + etc

\begin{cases} f_x(x,y,z) = z \\ g_x(x,y,z) = z + etc \end{cases}

K_{1y} = h f(x_i, y_i, z_i)

K_{1z} = hg(x_i, y_i, z_i)

K_{2y} = h f(x_i +h, y_i + K_{1y} , z_i + K_{1z})

K_{2z} = h g(x_i +h, y_i + K_{1y}, z_i + K_{1z})

y_{i+1}=y_i+\frac{K_{1y}+K_{2y}}{2}

z_{i+1}=z_i+\frac{K_{1z}+K_{2z}}{2}

x_{i+1} = x_i +h

import numpy as np

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras,
                   vertabla=False, precision = 6):
    ''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
        x0,y0 son valores iniciales, h es tamano de paso,
        muestras es la cantidad de puntos a calcular.
        f(x,y,z) = z #= y'
        g(x,y,z) = expresion con z=y'
    '''
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,3+4),dtype=float)

    # incluye el punto [x0,y0,z0,K1y,K1z,K2y,K2z]
    tabla[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    i=0
    if vertabla==True:
        print('Runge-Kutta Segunda derivada')
        print('i ','[ xi,  yi,  zi',']')
        print('   [ K1y,  K1z,  K2y,  K2z ]')
        np.set_printoptions(precision)
        print(i,tabla[i,0:3])
        print('  ',tabla[i,3:])
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
        if vertabla==True:
            txt = ' '
            if i>=10:
                txt = '  '
            print(str(i)+'',tabla[i,0:3])
            print(txt,tabla[i,3:])
    return(tabla)

[ Runge Kutta d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]
..

3. Runge-Kutta 4do Orden para Segunda derivada d²y/dx² en Python

import numpy as np

def rungekutta4_fg(fx,gx,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    tabla = np.zeros(shape=(tamano,3+8),dtype=float)

    # incluye el punto [x0,y0]
    tabla[0] = [x0,y0,z0,0,0,0,0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    
    for i in range(1,tamano,1):
        K1y = h * fx(xi,yi,zi)
        K1z = h * gx(xi,yi,zi)
        
        K2y = h * fx(xi+h/2, yi + K1y/2, zi + K1z/2)
        K2z = h * gx(xi+h/2, yi + K1y/2, zi + K1z/2)
        
        K3y = h * fx(xi+h/2, yi + K2y/2, zi + K2z/2)
        K3z = h * gx(xi+h/2, yi + K2y/2, zi + K2z/2)

        K4y = h * fx(xi+h, yi + K3y, zi + K3z)
        K4z = h * gx(xi+h, yi + K3y, zi + K3z)

        yi = yi + (K1y+2*K2y+2*K3y+K4y)/6
        zi = zi + (K1z+2*K2z+2*K3z+K4z)/6
        xi = xi + h
        
        tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z,K3y,K3z,K4y,K4z]
    return(tabla)

[ Runge Kutta d²y/dx² ] Algoritmo: [ RK 2do Orden ] [ RK 4to Orden] [ Ejercicio ]
..