Runge Kutta 4在C++中用於Lorenz系統

Question

有人可以在我的代碼中找到錯誤嗎？它可以在特殊的點上工作，但容差不適合該方法。我的錯誤是非常簡單的，所以我不認爲它的問題。請幫忙。

#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. 

作爲一個僅知道從計算RK4步驟的值y1=y+e和y2=y+e2=y+e/16，人們發現

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. 

這個公式確定的比循環快得多數千次迭代，每次有3個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