/*  prog4a.c -- Compute the Lorenz System
    Written Nov 12, 2002 by Eric Olson for Mathematics 483

    For information about the Lorenz System see section 4.3 in Applied 
    Analysis of the Navier Stokes Equations by Doering and Gibbon.

    For information about the Euler and RK methods used in solving this 
    equation see chapter 4 in Numerical Analysis for Applied Mathematics,
    Science and Engineering by Greenspan and Casulli.  */

#include <stdio.h>
#include <math.h>

#define N 3

double sigma=10.0;
double b=2.666666666666666666667;
double r=28.0;

void force(double F[N],double X[N]){
	F[0]=sigma*(X[1]-X[0]);
	F[1]=r*X[0]-X[1]-X[0]*X[2];
	F[2]=X[0]*X[1]-b*X[2];
	return;
}

void euler(double X[N],double tn,int n) {
	double F[3];
	double h=tn/n;
	int i,j;
	for(j=0;j<n;j++) {
		force(F,X);
		for(i=0;i<N;i++) {
			X[i]=X[i]+h*F[i];
		}
	}
	return;
}

void RK2(double X[N],double tn,int n){
	double Y[N];
	double K1[N],K2[N];
	double h=tn/n;
	int i,j;
	for(j=0;j<n;j++) {
		force(K1,X);
		for(i=0;i<N;i++) K1[i]*=h;
		for(i=0;i<N;i++) Y[i]=X[i]+K1[i];
		force(K2,Y);
		for(i=0;i<N;i++) K2[i]*=h;
		for(i=0;i<N;i++) X[i]=X[i]+0.5*(K1[i]+K2[i]);
	}
	return;
}

void RK4(double X[N],double tn,int n){
	double Y[N];
	double K1[N],K2[N],K3[N],K4[N];
	double h=tn/n;
	int i,j;
	for(j=0;j<n;j++) {
		force(K1,X);
		for(i=0;i<N;i++) K1[i]*=h; 
		for(i=0;i<N;i++) Y[i]=X[i]+0.5*K1[i];
		force(K2,Y);
		for(i=0;i<N;i++) K2[i]*=h;
		for(i=0;i<N;i++) Y[i]=X[i]+0.5*K2[i];
		force(K3,Y);
		for(i=0;i<N;i++) K3[i]*=h;
		for(i=0;i<N;i++) Y[i]=X[i]+K3[i];
		force(K4,Y);
		for(i=0;i<N;i++) K4[i]*=h;
		for(i=0;i<N;i++) X[i]=X[i]+(K1[i]+2*K2[i]+2*K3[i]+K4[i])/6.0;
	}
	return;
}

void init(double X[N]){
	X[0]=2.0; X[1]=3.0; X[2]=15.0; 
	return;
}

double norm(double X1[N],double X2[N]){
	double s=0.0;
	int i;
	for(i=0;i<N;i++) {
		double t=X1[i]-X2[i];
		s+=t*t;
	}
	return sqrt(s);
}

double converge(void (*solve)()){
	double tn=1.0;
	int n=64;
	double X[N],Xold[N];
	double en;
	int i,k;
	init(Xold);
	solve(Xold,tn,n);
	printf("n=%d, y(%lg)=[%lg, %lg, %lg]",
			n,tn,Xold[0],Xold[1],Xold[2]);
	for(k=0;k<10;k++){
		n*=2;
		init(X);
		solve(X,tn,n);
		en=norm(X,Xold);
		printf(", E=%lg\nn=%d, y(%lg)=[%lg, %lg, %lg]",
			en,n,tn,X[0],X[1],X[2]);
		for(i=0;i<N;i++) Xold[i]=X[i];
	}
	putchar('\n');
}

int main(){
	printf("prog4a.c -- Compute the Lorenz System\n");

	printf("\nEuler's Method:\n");
	converge(euler);

	printf("\nRK2 Method:\n");
	converge(RK2);

	printf("\nRK4 Method:\n");
	converge(RK4);

	return 0;
}