Runge-Kutta 4th order to solve 2nd order ODE using C++

Question

I tried to write a code to solve a 2nd order ODE but somehow it did not work as I intended.

the equation is 2y" + 5y' +3y = 0 ,y(0) = 3 and y'(0) = -4

The final answer will be y(x) = 3-4x << edit this is wrong

therefore y(1) = -1 and y'(1) = -4

But the program output is y(1) = 0.81414 y'(1) = -1.03727 << correct

Please help!

Thank you!

#include "stdafx.h"
#include <iostream>
#include <cmath>
#include <math.h> 

double ddx(double x,double y,double z)
{
    double dx = (-5*z-3*y)/2;
    return dx;
}

double test2ndorder(double x0, double y0, double z0, double xt, double h)
{
    int n = (int)((xt - x0) / h);
    double k1, k2, k3, k4;
    double l1, l2, l3, l4;

    double x = x0;
    double y = y0;
    double z = z0;

    for (int i = 1; i <= n; i++)
    {

        k1 = h * z;
        l1 = h * ddx(x, y, z);
        k2 = h * (z + 0.5*l1);
        l2 = h * ddx(x + 0.5*h, y + 0.5*k1, z + 0.5*l1);
        k3 = h * (z + 0.5*l2);
        l3 = h * ddx(x + 0.5*h, y + 0.5*k2, z + 0.5*l2);
        k4 = h * (z + l3);
        l4 = h * ddx(x + h, y + k3, z + l3);

        y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4);
        z = z + (1.0 / 6.0)*(l1 + 2 * l2 + 2 * l3 + l4);
        x = x + h;
        std::cout << y << "    ";
        std::cout << z << "\n";
    }
    return y;
}

int main()
{
    double x0, y0, z0, x, y, z,h;
    x0 = 0;
    x = 1;
    y0 = 3;
    z0 = -4;
    h =0.01;

    y = test2ndorder(x0, y0, z0, x, h);
    std::cout << y;
}

Answer 1

Code in C++ language. In this implementation the following initial conditions were used, x0 = 0, y(x0) = 3 and u(x0) = -4.

// Autor    : Carlos Eduardo da Silva Lima
// Tema     : Runge-Kutta de quarta ordem
// Linguagem: C++ (ANSI)
// IDE      : Online GBD
// Data     : 19/11/2022

#include <iostream>
#include <cmath>

#define N 10000

float edo1(float x, float y, float u);
float edo2(float x, float y, float u);
void RK4(float x0, float y0, float u0, float h);

int main()
{
    
    RK4(0, 3, -4, 1E-3);

    return 0;
}


float edo1(float x, float y, float u){
    return u;
}

float edo2(float x, float y, float u){
    return ((-3*y-5*u)/2);
}

void RK4(float x0, float y0, float u0, float h){
    
    float x[N], y[N], u[N];
    
    x[0] = x0;
    y[0] = y0;
    u[0] = u0;
    
    int i = 0;
    do{
        
        float k11 = h*edo1(x[i],y[i],u[i]);
        float k12 = h*edo2(x[i],y[i],u[i]);
        
        float k21 = h*edo1(x[i]+(h/2),y[i]+(k11/2),u[i]+(k12/2));
        float k22 = h*edo2(x[i]+(h/2),y[i]+(k11/2),u[i]+(k12/2));
        
        float k31 = h*edo1(x[i]+(h/2),y[i]+(k21/2),u[i]+(k22/2));
        float k32 = h*edo2(x[i]+(h/2),y[i]+(k21/2),u[i]+(k22/2));
        
        float k41 = h*edo1(x[i]+h,y[i]+k31,u[i]+k32);
        float k42 = h*edo2(x[i]+h,y[i]+k31,u[i]+k32);
        
        y[i+1] = y[i] + ((k11+2*(k21+k31)+k41)/6);
        u[i+1] = u[i] + ((k12+2*(k22+k32)+k42)/6);
        x[i+1] = x[i] + h;
        
        i++;
    }while(i<=(N-1));
    
    // Saída
    int j = 0;
    do{
        std::cout << "x = " << x[j] << " | y = " << y[j] << " | u = " << u[j] << std::endl;
        j++;
    }while(j<=(N-1));
}

I also did the implementation in the python language. Here I use an integration method for edo's, developed and documented in scipy, scipy.integrate.odeint. If you don't agree, or if you find any logic errors, please let me know. Up until :).

# Autor    : Carlos Eduardo da Silva Lima
# Tema     : Runge-Kutta de quarta ordem
# Linguagem: Python
# IDE      : Google Colab
# Data     : 19/11/2022

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Condições iniciais 
x_inicial = 0.0
x_final = 10.0
y0 = 3
u0 = -4
N = 10000

# Equações diferenciais na primeira ordem
def edo(r,x):
  y = r[0]
  u = r[1]
  return np.array([u, (-3*y-5*u)/2])

# Aplicação do métodos para resolução das edo´s declaradas acima emonjunto com as condições
x = np.linspace(x_inicial, x_final, N)
r0 = np.array([y0,u0])

# Resolução 
sol = odeint(edo, r0, x)

# Grvandos os resultados em novas variáveis
y = sol[:,0]
u = sol[:,1]

# Plot
plt.style.use('dark_background')
plt.figure(figsize=(7,7))
plt.plot(x,y,'m-', linewidth = 3.5)
plt.grid(color='y', linestyle='-', linewidth=0.05)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Odeint python")
plt.show()

Answer 2

2y" + 5y' +3y = 0 ,y(0) = 3 and y'(0) = -4

has characteristic roots that are solutions of

(2r)^2 + 5*(2r) + 6 = 0  <==>  r = -1 or r = -1.5

so that the solution is

y(x) = A*exp(-x)+B*exp(-1.5*x)
 3 = y(0)  =  A +     B
-4 = y'(0) = -A - 1.5*B

so that B = 2 and A = 1 which gives

y(1)  =  0.814139761468302 
y'(1) = -1.03726992161673

which is about what you got.

Answer 3

In your for loop in test2ndorder, you use the constant x0, y0, and z0 values, so the values computed in each iteration won't change. You also compute a new value for y and never use it.

It appears that you should be using x, y, and z instead.