численные методы — Рунге Кутта 4 в C ++ для системы Лоренца

Может кто-нибудь найти ошибку в моем коде? Это работает нормально в специальных точках, но допуск не подходит для метода. Мой 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);
}

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

Эта формула намного быстрее, чем цикл из нескольких тысяч итераций с 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
1

Другие решения

Других решений пока нет …

По вопросам рекламы ammmcru@yandex.ru
Adblock
detector