Solving and plotting the Lorenz Attractor

plot.py

#!/usr/bin/env python3

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt

if __name__ == "__main__":
    with open("lorenz.txt","r",encoding='utf-8') as f:
        tmp=[i.split("   ") for i in f.readlines()]
    
    x=[float(i[0]) for i in tmp]
    y=[float(i[1]) for i in tmp]
    z=[float(i[2]) for i in tmp]
    
    xs=np.array(x)
    ys=np.array(y)
    zs=np.array(z)
    
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot(xs, ys, zs)
    
    plt.savefig("lorenz.png", format="png")

runge_kutta.c

#include "runge_kutta.h"

double * rk4(double *z0, double tau, double * (*dgl)(), double T, int dim)
{
    double *zout, *z_seq, t;
    int j, n, N;

    N = T/tau;

    z_seq = (double *) calloc(dim * N, sizeof(double));
    zout = (double *) calloc(dim, sizeof(double));

    for(j=0;j<dim;j++)
        z_seq[j * N + 0] = z0[j];

    for(n=1;n<N;n++)
    {
        t=n*tau;
        for(j=0;j<dim;j++)
            zout[j] = z_seq[j * N + n-1];

        zout = rk4_step(zout,t,tau,dgl,dim);

        for(j=0;j<dim;j++)
            z_seq[j * N + n] = zout[j];
    }
    return z_seq;
}

double * rk4_step(double *z, double t, double tau, double * (*dgl)(), int dim)
{
    int i;
    double *ztmp, *Hs[4];

    for(i=0;i<4;i++)
    {
        Hs[i] = (double *) calloc(dim, sizeof(double));
    }
    ztmp = (double *) calloc(dim, sizeof(double));

    Hs[0]=dgl(z, t, dim);
    for(i=0;i<dim;i++)
        ztmp[i] = z[i]+tau/2*Hs[0][i];
    Hs[1]=dgl(ztmp, t+tau/2, dim);
    for(i=0;i<dim;i++)
        ztmp[i] = z[i]+tau/2*Hs[1][i];
    Hs[2]=dgl(ztmp, t+tau/2, dim);
    for(i=0;i<dim;i++)
        ztmp[i] = z[i]+tau  *Hs[2][i];
    Hs[3]=dgl(ztmp, t+tau, dim);

    for(i=0;i<dim;i++)
        ztmp[i] = z[i] + tau/6*(Hs[0][i]+2*Hs[1][i]+2*Hs[2][i]+Hs[3][i]);

    return ztmp;
}

double * rk_lorenz_func(double *r, double t, int dim)
{
    // Lorenz Modell
    double sigma, r_param, b;
    double x, y, z, *out;
    out = (double *) calloc(dim, sizeof(double));
    sigma   = 10;
    r_param = 28;
    b       = 8/3;

    x=r[0]; y=r[1]; z=r[2];

    out[0] = sigma*(y-x);
    out[1] = r_param*x - y - x*z;
    out[2] = x*y - b*z;

    return out;
}


void runge_kutta_lorenz()
{
    double z0[]={1,1,20}, tau = 0.01, T = 20, *z_seq;
    int dim=3, j, i, N;

    N = T/tau;

    z_seq = rk4(z0, tau, rk_lorenz_func, T, dim);
    for(i=0;i<N;i++)
    {
        for(j=0;j<dim;j++)
            printf("% 3.6g   ",z_seq[j * N + i]);
        printf("\n");
    }
}

void main()
{
    runge_kutta_lorenz();
}

runge_kutta.h

#ifndef RUNGE_KUTTA_H
#define RUNGE_KUTTA_H

#include <stdio.h>
#include <stdlib.h>

double * rk4_step(double *z, double t, double tau, double * (*dgl)(), int dim);
double * H(double *z, double t, int dim);

double * rk_lorenz_func(double *r, double t, int dim);

#endif //RUNGE_KUTTA_H