Skip to content

Instantly share code, notes, and snippets.

@surt91
Created Sep 16, 2016
Embed
What would you like to do?
Solving and plotting the Lorenz Attractor
#!/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")
#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();
}
#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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment