[Python] Why my c++ code is slower than python code doing same thing?

I'm new to the programming. I've tried to make 3 scripts in Python, c++ and Matlab that would calculate Voltage over Inductor in RL circuit using matrices and finite difference method. When i ran those codes it appeared that MatLab took 0.4s to compute, Python took 1s to compute and c++ over 30s. Here I'm posting c++ code:

#include <iostream>
#include <vector>
#include <cmath>
#include <chrono>
#include <Eigen/Dense>
#include "matplotlibcpp.h"
#include <omp.h>

namespace plt = matplotlibcpp;

double R = 10.0;
double L = 2.0;
double V0 = 5.0;
double t1 = 1.0;
double t2 = 2.0;
double t_start = 0.0;
double t_end = 3.0;

double exact_voltage(double t) {
if (t < t1) {
return 0;
} else if (t1 <= t && t < t2) {
return V0 * std::exp(-(R / L) * (t - t1));
} else {
return -V0 * std::exp(-(R / L) * (t - t2));
}
}

void solve(double dt, std::vector<double>& solved_s_times, std::vector<double>& solved_V_L) {
int N = static_cast<int>((t_end - t_start) / dt);

Eigen::MatrixXd A = Eigen::MatrixXd::Zero(N, N);
Eigen::VectorXd B = Eigen::VectorXd::Zero(N);

for (int i = 0; i < N; ++i) {
A(i, i) = L / dt + R;
if (i > 0) {
A(i, i - 1) = -L / dt;
}
}

std::vector<double> solve_times(N);
for (int i = 0; i < N; ++i) {
solve_times = t_start + i * dt;
if (t1 <= solve_times && solve_times < t2) {
B(i) = V0;
} else if (solve_times >= t2) {
B(i) = 0.0;
}
}

Eigen::VectorXd I = A.colPivHouseholderQr().solve(B);

Eigen::MatrixXd A_dI_dt = Eigen::MatrixXd::Zero(N, N);
for (int i = 0; i < N - 1; ++i) {
A_dI_dt(i, i) = -1.0 / dt;
A_dI_dt(i, i + 1) = 1.0 / dt;
}

Eigen::VectorXd dI_dt = A_dI_dt * I;

solved_V_L.resize(N);
solved_s_times.resize(N);

for (int i = 0; i < N; ++i) {
solved_V_L = L * dI_dt;
solved_s_times = solve_times + dt;
}
}

int main() {

std::cout << Eigen::nbThreads() << std::endl;
std::vector<double> dt_values = {0.002, 0.001, 0.004};

auto start = std::chrono::high_resolution_clock::now(); // Start pomiaru czasu

for (double dt : dt_values) {
std::vector<double> times, V_L;
solve(dt, times, V_L);
plt:lot(times, V_L, {{"label", "Symulacja (dt = " + std::to_string(dt) + ")"}});
}

auto end = std::chrono::high_resolution_clock::now(); // Koniec pomiaru czasu
std::chrono::duration<double> duration = end - start;
std::cout << "Wykonano w: " << duration.count() << "s" << std::endl;

std::vector<double> exact_times, exact_V_L;
double dt_exact = 0.001;
int N_exact = static_cast<int>((t_end - t_start) / dt_exact);

exact_times.resize(N_exact);
exact_V_L.resize(N_exact);

for (int i = 0; i < N_exact; ++i) {
exact_times = t_start + i * dt_exact;
exact_V_L = exact_voltage(exact_times);
}

plt:lot(exact_times, exact_V_L, {{"label", "Dokładne rozwiązanie"}});
plt::xlabel("Czas (s)");
plt::ylabel("Napięcie (V)");
plt::title("Napięcie na cewce w obwodzie RL - symulacja i dokładne rozwiązanie");
plt::grid(true);
plt::legend();
plt::show();

return 0;
}


And the Python code:

import numpy as np
import matplotlib.pyplot as plt
import time
from numba import njit

R = 10.0
L = 2.0
V0 = 5.0

t_start = 0.0
t_end = 3
dt_values = [0.002, 0.001, 0.004]

plt.figure(figsize=(12.8, 7.2))

t1 = 1
t2 = 2


# @njit
def exact_voltage(t):
if t < t1:
return 0
elif t1 <= t < t2:
return V0 * np.exp(-(R / L) * (t - t1))
else:
return -V0 * np.exp(-(R / L) * (t - t2))


# @njit
def solve(dt):
solve_times = np.arange(t_start, t_end, dt)
N = len(solve_times)

A = np.zeros((N, N))
B = np.zeros(N)

for i in range(N):
A[i, i] = L / dt + R
if i > 0:
A[i, i - 1] = -L / dt

for i in range(N):
if t1 <= solve_times < t2:
B = V0
elif solve_times >= t2:
B = 0.0

I = np.linalg.solve(A, B)

A_dI_dt = np.zeros((N, N))

for i in range(N - 1):
A_dI_dt[i, i] = - 1 / dt
A_dI_dt[i, i + 1] = 1 / dt

dI_dt = np.dot(A_dI_dt, I)

solved_V_L = L * dI_dt

solved_s_times = np.zeros(N)
for i in range(N):
solved_s_times = solve_times + dt

return solved_s_times, solved_V_L


start = time.time()
for dt_value in dt_values:
times, V_L = solve(dt_value)
plt.plot(times, V_L, label=f'Symulacja (dt = {dt_value})')
end = time.time()
print(f"Wykonano w: {end - start}s")
exact_times = np.arange(t_start, t_end, 0.001)
exact_V_L = [exact_voltage(t) for t in exact_times]
plt.plot(exact_times, exact_V_L, 'k--', label='Dokładne rozwiązanie')

plt.xlabel('Czas (s)')
plt.ylabel('Napięcie (V)')
plt.title('Napięcie na cewce w obwodzie RL - symulacja i dokładne rozwiązanie')
plt.grid(True)
plt.legend()
plt.show()


I've tried using different c++ compilers, adding -O3 and -fopenmp to build. I added & to functions too when i watched video on youtube explaining poor c++ performance. None of these actions helped much, they all subtracted around 1s from computation time of c++ code. Maybe someone sees something obvious that i'm missing? I'll be very thankful for explaining why my code is that much slower than Python and MatLab. If you could suggest any changes to the code that will help speed it up it would be great.

Continue reading...