Метод Рунге-Кутта 4-го порядка

@Luckyscissors · Регистрация: 30.10.2016

Правильно ли я записал выражения для коэффициентов моего уравнения(система Лоренца)

C++

#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 deltapodshet(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 prirash;
    
    prirash.dx=a;
    prirash.dy=b;
    prirash.dz=c;
    
    return prirash;
    
}
int main ()
{
 
 
    double dx,dy,dz,h,ei,xi,zi,yi,e,t=0,T=0.08,a,b,c,k=0,E;
    int n=0,m=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<<"Pogreshost=";
    cin>>E;
    xi=x;
    
do
    {
    
    
    delta konec1=deltapodshet(xi,yi,zi,h/2);
    a=konec1.dx;
    b=konec1.dy;
    c=konec1.dz;
    konec1=deltapodshet(a,b,c,h/2);
    e=fabs((x-konec1.dx)/15);
 
    
    while (e>E)
    {
    h=h/2;
    delta konec3=deltapodshet(x,y,z,h);
    xi=konec3.dx;
    konec3=deltapodshet(x,y,z,h/2);
    a=konec3.dx;
    b=konec3.dy;
    c=konec3.dz;
    konec3=deltapodshet(a,b,c,h/2);
    e=fabs((xi-konec3.dx)/15); 
    };
    xi=x;
    yi=y;
    zi=z;
    ei=e;
    
    konec1=deltapodshet(x,y,z,h);
    
    x=konec1.dx;
    y=konec1.dy;
    z=konec1.dz;
    
    t+=h;   
    k+=1;
//  cout<<"x="<<x<<" y= "<<y<<" z="<<z<<" Pogreshnost= "<<e<<" Time="<<t<<endl;
    fprintf(fff,"%lf, %lf, %lf, %lf ,%lf ,%lf\n", x, y,z,e,t, h);               
    
}
while (t<=T);
fprintf(fff,"%d",k); 
fclose(fff);
}

Кажется , я все делал правильно , но моя погрешность постоянно на несколько порядков меньше нужной ( более того она колеблется :может увеличиваться и уменьшаться)

@afront · 10.12.2016, 21:01

C++

#include <iostream>
#include <cmath>
 
using namespace std;
 
static double f(double x, double y)
{
    return y*cos(x);
}
double ODE_RungeKutta4(double x0, double y0, double h, double x)
{
    double xnew, ynew, k1, k2, k3, k4, result = 0;
    if (x == x0)
        result = y0;
    else if (x > x0)
    {
        do
        {
            if (h > x - x0) h = x - x0;
            k1 = h * f(x0, y0);
            k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1);
            k3 = h * f(x0 + 0.5 * h, y0 + 0.5 * k2);
            k4 = h * f(x0 + h, y0 + k3);
            ynew = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            xnew = x0 + h;
            x0 = xnew;
            y0 = ynew;                                        
        }while (x0 < x);
        result = ynew;
    }
    return result;
}
int main ()
{
    double h = 0.001;
    double x0 = 0.0;
    double y0 = 1.0;
    double result = y0;
    for (int i = 0; i < 11; i++)
    {
        double x = 0.1 * i;
        result = ODE_RungeKutta4(x0, result, h, x);
        double exact = exp(sin(x));
        if (i % 5 == 0)
            cout<<" x = "<<x<<" y = "<< result<<" exact = "<<exact<<endl;
        x0 = x;
    }
    system("pause");
    return 0;
}

@Luckyscissors · 10.12.2016, 21:29 **[ТС]**

Спасибо , но это немного не то , что мне нужно. Если бы у вас нашелся пример метода Рунге-Кутта для системы из трех дифференциальных уравнений , я был бы очень благодарен.

@afront · 11.12.2016, 09:57

Luckyscissors, у меня есть готовый пример решения системы уравнений осциллятора Лоренца методом Эйлера

C++

#include <iostream>
 
using namespace std;
 
static double dx(double x, double y, double z)
{ return 10.0 * (y - x); }
 
static double dy(double x, double y, double z)
{ return x * (28.0 - z) - y; }
 
static double dz(double x, double y, double z)
{ return x * y - 8.0 * z / 3.0; }
 
int main()
{
    double x = 0.0, y = 5.0, z = 10.0;  //initial conditions
    double dt = 0.001;                  //step size
 
    for (int i = 0; i < 100; i++)
    {
        double xnew = x + dx(x, y, z) * dt;
        double ynew = y + dy(x, y, z) * dt;
        double znew = z + dz(x, y, z) * dt;
        x = xnew;
        y = ynew;
        z = znew;
        cout<<" x = "<<x<<" y = "<< y<<" z = "<< z<<endl;
    }
system("pause");
return 0;
}

@Luckyscissors · 11.12.2016, 12:38 **[ТС]**

Спасибо. А не могли бы вы еще подсказать принцип работы автовыбора шага ? Если вы ,конечно, не против. У меня получилось реализовать только "автовыбор по убыванию"(не знаю как назвать), он будет только уменьшать шаг и если я введу уже изначально маленький шаг ,то автовыбор его не изменит и программа будет считать с ним.

@afront · 11.12.2016, 14:48

посмотрите эту статью, там есть что то по этому поводу
https://www.codeproject.com/Ar... -equations
а в этой книжке расписан метод Мерсона
http://kachat-knigi.ru/excel-u... -VBA-C.htm

@Luckyscissors · 12.12.2016, 14:11 **[ТС]**

Извините , за некоторую назойливость , но не могли бы вы проверить мою программу ? Я прикрутил туда довольно топорный автовыбор шага ( на основе текста учебника ,еще я прочитал статью там использовался try_stepper, который я не очень понимаю , как реализовать , но его смысл понятен). Графики в особых точках соответствуют теории , но с погрешностями все равно беда.

+табличка для метода автовыбора шага не соответствует действительности ( c уменьшением порядка погрешности в 4 раза шаг должен уменьшаться на один порядок)

C++

#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 deltapodshet(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 prirash;
    
    prirash.dx=a;
    prirash.dy=b;
    prirash.dz=c;
    
    return prirash;
    
}
double Error(double x,double y ,double z, double h)
{
    double xi,e,a,b,c;
    delta konec3=deltapodshet(x,y,z,h);
    xi=konec3.dx;
    konec3=deltapodshet(x,y,z,h/2);
    a=konec3.dx;
    b=konec3.dy;
    c=konec3.dz;
    konec3=deltapodshet(a,b,c,h/2);
    e=fabs((xi-konec3.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<<"Pogreshost=";
    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<<"vse 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=deltapodshet(x,y,z,h);
    
    x=konec1.dx;
    y=konec1.dy;
    z=konec1.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,"Средний шаг = %lf",T/k); 
fclose(fff);
}

@afront · 12.12.2016, 15:41

Luckyscissors, выше я вам дал решение примерно той же задачи но другим методом, сравните с ней, подставте оттуда числа в свой код и сравните результаты

@Luckyscissors · 12.12.2016, 23:15 **[ТС]**

На заданном в вашем решении промежутке времени различий между моим методом и вашим нету. Но как я понимаю , предоставленное решение методом Эйлера , сделано с целью увидеть аттрактор Лоренца ? В моей программе при r=28(когда
должен наблюдаться аттрактор) выходит какая-то мишура из точек (строю в Origin XYZ график для фазового портрета) , причем x(t), y(t) притягиваются к сложному множеству аттрактору.

При промежутке времени в 100 секунд , у меня, как я уже сказал мишура , а в методе Эйлера хорошо виден аттрактор.

Добавлено через 52 минуты
Я выкинул все ,что связано с вычислением погрешности и оставил только счет ( как в методе Эйлера), и все заработало.
Чувствую некоторое облегчение , что хотя бы коэф я записал правильно.

C++

 double x = 0.0, y = 5.0, z = 10.0;  //initial conditions
    double dt = 0.001;                  //step size
 
    for (int n = 0; n < 10000; n++)
    {
        double k1,k2,k3,k4,m1,m2,m3,m4,p1,p2,p3,p4,dx,dy,dz,a,b,c;
    
    
    k1=dt*f1(x,y);
    m1=dt*f2(x,y,z);
    p1=dt*f3(x,y,z);
        
    
    k2=dt*f1((x+k1/2),(y+m1/2));
    m2=dt*f2((x+k1/2) ,(y+m1/2),(z+p1/2));
    p2=dt*f3((x+k1/2),(y+m1/2),(z+p1/2));
    
    
    k3=dt*f1((x+k2/2),(y+m2/2));
    m3=dt*f2((x+k2/2),(y+m2/2),(z+p2/2));
    p3=dt*f3((x+k2/2),(y+m2/2),(z+p2/2));
    
 
    k4=dt*f1((x+k3),(y+m3));
    m4=dt*f2((x+k3),(y+m3),(z+p3));
    p4=dt*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;
    
    x+=dx;
    y+=dy;
    z+=dz;
    
    fprintf(fff,"%lf %lf %lf\n",x,y,z);
       // cout<<" x = "<<x<<" y = "<< y<<" z = "<< z<<endl;
       k=k+1;
    }
 
 
fprintf(fff,"Step = %lf",100/k); 
fclose(fff);
}

@afront · 13.12.2016, 07:55

можно продолжить сравнения, подставить теперь в метод Эйлера ваши данные и посмотреть что получится

@Luckyscissors · 17.12.2016, 17:54 **[ТС]**

На больших промежутках времени даже предоставленное решение методом Эйлера теряет форму аттрактора Лоренца ( который по сути должен быть усточивым и принимать определенную форму при t->беск.). Может ли это быть ограничением метода ?

@afront · 17.12.2016, 18:03

метод Эйлера первого порядка, может быть такое, мне кажется

@Luckyscissors · 18.12.2016, 23:07 **[ТС]**

Извините , тут обнаружилась еще одна ошибка. В готовом решении и моем , при r=19 наблюдаются неустойчивые фокусы(аттрактор Лоренца) , хотя при данном r должны получаться устойчивые фокусы.Аттрактор появляется при пороговом r>24,7. Исключая ошибки в теории , остаются только ошибки в коде. Вы можете предположить , что может быть не так ?

@afront · 19.12.2016, 09:14

попробуйте проверить это в еще какой то реализации метода, например в этой
https://www.codeproject.com/Ar... -equations

@Luckyscissors 0 / 0 / 0 Регистрация: 30.10.2016 Сообщений: 8
	10.12.2016, 21:29 [ТС]	3
	Спасибо , но это немного не то , что мне нужно. Если бы у вас нашелся пример метода Рунге-Кутта для системы из трех дифференциальных уравнений , я был бы очень благодарен. 0

@Luckyscissors 0 / 0 / 0 Регистрация: 30.10.2016 Сообщений: 8
	11.12.2016, 12:38 [ТС]	5
	Спасибо. А не могли бы вы еще подсказать принцип работы автовыбора шага ? Если вы ,конечно, не против. У меня получилось реализовать только "автовыбор по убыванию"(не знаю как назвать), он будет только уменьшать шаг и если я введу уже изначально маленький шаг ,то автовыбор его не изменит и программа будет считать с ним. 0

@afront 1481 / 1198 / 819 Регистрация: 29.02.2016 Сообщений: 3,579
	11.12.2016, 14:48	6
	посмотрите эту статью, там есть что то по этому поводу https://www.codeproject.com/Ar... -equations а в этой книжке расписан метод Мерсона http://kachat-knigi.ru/excel-u... -VBA-C.htm 1

@afront 1481 / 1198 / 819 Регистрация: 29.02.2016 Сообщений: 3,579
	12.12.2016, 15:41	8
	Luckyscissors, выше я вам дал решение примерно той же задачи но другим методом, сравните с ней, подставте оттуда числа в свой код и сравните результаты 1

@afront 1481 / 1198 / 819 Регистрация: 29.02.2016 Сообщений: 3,579
	13.12.2016, 07:55	10
	можно продолжить сравнения, подставить теперь в метод Эйлера ваши данные и посмотреть что получится 1

@Luckyscissors 0 / 0 / 0 Регистрация: 30.10.2016 Сообщений: 8
	17.12.2016, 17:54 [ТС]	11
	На больших промежутках времени даже предоставленное решение методом Эйлера теряет форму аттрактора Лоренца ( который по сути должен быть усточивым и принимать определенную форму при t->беск.). Может ли это быть ограничением метода ? 0

@afront 1481 / 1198 / 819 Регистрация: 29.02.2016 Сообщений: 3,579
	17.12.2016, 18:03	12
	метод Эйлера первого порядка, может быть такое, мне кажется 1

@Luckyscissors 0 / 0 / 0 Регистрация: 30.10.2016 Сообщений: 8
	18.12.2016, 23:07 [ТС]	13
	Извините , тут обнаружилась еще одна ошибка. В готовом решении и моем , при r=19 наблюдаются неустойчивые фокусы(аттрактор Лоренца) , хотя при данном r должны получаться устойчивые фокусы.Аттрактор появляется при пороговом r>24,7. Исключая ошибки в теории , остаются только ошибки в коде. Вы можете предположить , что может быть не так ? 0

@afront 1481 / 1198 / 819 Регистрация: 29.02.2016 Сообщений: 3,579
	19.12.2016, 09:14	14
	попробуйте проверить это в еще какой то реализации метода, например в этой https://www.codeproject.com/Ar... -equations 1

Опции темы

Метод Рунге-Кутта 4-го порядка

Решение