Can someone find the error in my code ? It work ok in the special dots , but the tolerance doesn't fit to the method. My Errorstepper is pretty simple , so i don't think that the problem in it. Please help.
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <time.h>
#include <iostream>
#include <string>
using namespace std;
const int q=10;
const double b=8/3;
double r;
double x,y,z;
struct delta
{
double dx;
double dy;
double dz;
};
double f1(double x , double y)
{
return (q*(y-x));
}
double f2(double x,double y, double z)
{
return ((-x)*z+r*x-y);
}
double f3(double x,double y,double z)
{
return (x*y- b*z);
}
delta calc(double x,double y,double z,double dh=0.1)
{
double k1,k2,k3,k4,m1,m2,m3,m4,p1,p2,p3,p4,dx,dy,dz,a,b,c,h;
h=dh;
k1=h*f1(x,y);
m1=h*f2(x,y,z);
p1=h*f3(x,y,z);
//h+=dh/2;
k2=h*f1((x+k1/2.),(y+m1/2.));
m2=h*f2((x+k1/2.) ,(y+m1/2.),(z+p1/2.));
p2=h*f3((x+k1/2.),(y+m1/2.),(z+p1/2.));
// h+=dh/2;
k3=h*f1((x+k2/2.),(y+m2/2.));
m3=h*f2((x+k2/2.),(y+m2/2.),(z+p2/2.));
p3=h*f3((x+k2/2.),(y+m2/2.),(z+p2/2.));
// h+=dh;
k4=h*f1((x+k3),(y+m3));
m4=h*f2((x+k3),(y+m3),(z+p3));
p4=h*f3((x+k3),(y+m3),(z+p3));
dx=(k1+2*k2+2*k3+k4)/6.;
dy=(m1+2*m2+2*m3+m4)/6.;
dz=(p1+2*p2+2*p3+p4)/6.;
a=x+dx;
b=y+dy;
c=z+dz;
delta add;
add.dx=a;
add.dy=b;
add.dz=c;
return add;
}
double Error(double x,double y ,double z, double h)
{
double xi,e,a,b,c;
delta end1=calc(x,y,z,h);
xi=end1.dx;
end1=calc(x,y,z,h/2);
a=end1.dx;
b=end1.dy;
c=end1.dz;
end1=calc(a,b,c,h/2);
e=fabs((xi-end1.dx)/15);
return e;
}
int main ()
{
double dx,dy,dz,h,ei,xi,zi,yi,e,t=0,T=0.08,a,b,c,k=0,E,h1,e1,m,E1;
int n=0;
FILE *fff;
fff=fopen("data.txt","w");
cout <<"x0=";
cin>>x;
cout<<"y0=";
cin>>y;
cout<<"z0=";
cin>>z;
cout<<"h=";
cin>>h;
cout<<"r=";
cin>>r;
cout<<"Time=";
cin>>T;
cout<<"Toleranse=";
cin>>E;
xi=x;
double hmin=0.00005;
if (E<=0.0000001)
hmin=hmin/100;
else cout<<"ok"<<endl;
E1=E/10;
do
{
e=Error(x,y,z,h);
while (e<E1 || e>E)
{
if (e<E1)
h+=hmin;
else if (e>E)
h-=hmin;
else cout<<" ok"<<endl;
e=Error(x,y,z,h);
};/*
while(e>E)
{
h=h/2;
e=Error(x,y,z,h);
};*/
xi=x;
yi=y;
zi=z;
ei=e;
delta konec1=calc(x,y,z,h);
x=end1.dx;
y=end1.dy;
z=end1.dz;
t+=h;
k+=1;
// cout<<"x="<<x<<" y= "<<y<<" z="<<z<<" Pogreshnost= "<<e<<" Time="<<t<<endl;
fprintf(fff,"%lf, %lf, %lf, %.10lf ,%lf ,%lf\n", x, y,z,e,t, h);
}
while (t<=T);
fprintf(fff,"Step = %lf",T/k);
fclose(fff);
}
Your error stepper, as simple as it seems, is implementing some wrong ideas.
The general assumption is that in lowest order the local error is e=C·h^5. Doing two steps with half the step size thus gives an error
e2=2·C·(h/2)^5=C·h^5/16=e/16.
As one only knows the values y1=y+e and y2=y+e2=y+e/16 from computing the RK4 steps, one finds that
y1-y2=15/16*e or e=16/15*(y1-y2) and e2=(y1-y2)/15.
Assuming that the local error uniformly and additively translates to the global error, one gets a global error of e·(T/h). Or interpreted differently, e/h is the local error density, the desired global error E translates into an average error density E/T. Thus you have to compare e to E/T·h, not to a naked E. Especially the missing factor h will result in unsuitable step sizes.
The whole step size computation also is computationally too expensive. As the local error has been found as e=C·h^5, and the desired local error is E·h/T, then the probably best step size h1 (excluding large higher order effects) is determined by
C·h1^4=E/T or h1 = h*(E/(e*T))^0.25.
This formula is much faster than a loop of several thousand iterations with 3 RK4 steps each. One can build checks around that formula so that one is sure that the new step size really has that local error and that the change in step size is not too radical, in pseudo-concept-C code
fac = 1.0;
do {
h *= fac;
y1 = /* RK4 step with h */;
y2 = /* 2 RK4 steps with h/2 */;
e = 16*norm(y1-y2)/15;
fac = pow((E*h)/(e*T), 0.25);
fac = min(20, max(0.05, fac) );
} while( fac<0.5 or 2<fac )
// advance the integration