Runge Kutta 4 в C++ для системы Lorenz

Question

Может кто-нибудь найти ошибку в моем коде? Он работает нормально в специальных точках, но допуск не подходит для метода. Мой ErrorStepper довольно прост, поэтому я не думаю, что проблема в нем. Пожалуйста помоги.

#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); 
} 

Answer 1

Ваш шатер ошибки, как ни странно, реализует некоторые неправильные идеи.

Общее предположение состоит в том, что в нижнем порядке локальная ошибка e=C·h^5. Делая два шага с половиной шага, таким образом, дает ошибку

e2=2·C·(h/2)^5=C·h^5/16=e/16. 

Как один знает только значения y1=y+e и y2=y+e2=y+e/16 из расчета шагов RK4, можно обнаружить, что

y1-y2=15/16*e or e=16/15*(y1-y2) and e2=(y1-y2)/15. 

В предположении, что локальная ошибка равномерно и аддитивно переводит на глобальную ошибку, появляется глобальная ошибка e·(T/h). Или интерпретируется иначе, e/h - локальная плотность ошибок, требуемая глобальная ошибка E означает среднюю плотность ошибок E/T. Таким образом, вы должны сравнить e с E/T·h, а не с голым E. Особенно недостающий фактор h приведет к неподходящим размерам шага.

---

Весь расчет размера шага также является слишком дорогостоящим. По мере того как локальная ошибка была найдена как e=C·h^5, и желаемый локальная ошибка E·h/T, то, вероятно, лучше размер шага h1 (за исключением больших высших эффектов порядка) определяется

C·h1^4=E/T or h1 = h*(E/(e*T))^0.25. 

Эта формула значительно быстрее, чем цикл из нескольких тысяч итераций с тремя шагами RK4 каждый. Можно построить проверки вокруг этой формулы, так что один уверен, что новый размер шага действительно имеет такую ​​локальную ошибку и что изменение размера шага не слишком радикально, в коде псевдо-концепция-C

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/(e*T), 0.25); 
    fac = min(20, max(0.05, fac)); 
} while(fac<0.5 or 2<fac) 
// advance the integration