0 / 0 / 0
Регистрация: 30.10.2016
Сообщений: 8
1

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

10.12.2016, 19:20. Показов 6230. Ответов 13
Метки нет (Все метки)

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

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#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);
}
Кажется , я все делал правильно , но моя погрешность постоянно на несколько порядков меньше нужной ( более того она колеблется :может увеличиваться и уменьшаться)
__________________
Помощь в написании контрольных, курсовых и дипломных работ, диссертаций здесь
0
Лучшие ответы (1)
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
10.12.2016, 19:20
Ответы с готовыми решениями:

Метод Рунге-Кутта 4 порядка
Помогите найти ошибку в методе рунге-Кутта 4 порядка System::System(double m, const Vector3D&amp; g,...

Метод Рунге-Кутта 2-го порядка
Написал код программы. С компиляцией вроде бы нет проблем. А цикл for воспринимать не хочет....

Метод Рунге-Кутта 3 порядка
Начерикал что-то,вроде бы работает,но не уверен,да и не доходит,как под условия сделать...

Метод Рунге-Кутта 4-го порядка
#include &lt;stdio.h&gt; #include &lt;iostream&gt; #include &lt;math.h&gt; using namespace std; double f(double...

13
1481 / 1198 / 819
Регистрация: 29.02.2016
Сообщений: 3,579
10.12.2016, 21:01 2
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#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;
}
2
0 / 0 / 0
Регистрация: 30.10.2016
Сообщений: 8
10.12.2016, 21:29  [ТС] 3
Спасибо , но это немного не то , что мне нужно. Если бы у вас нашелся пример метода Рунге-Кутта для системы из трех дифференциальных уравнений , я был бы очень благодарен.
0
1481 / 1198 / 819
Регистрация: 29.02.2016
Сообщений: 3,579
11.12.2016, 09:57 4
Лучший ответ Сообщение было отмечено Luckyscissors как решение

Решение

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

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#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;
}
2
0 / 0 / 0
Регистрация: 30.10.2016
Сообщений: 8
11.12.2016, 12:38  [ТС] 5
Спасибо. А не могли бы вы еще подсказать принцип работы автовыбора шага ? Если вы ,конечно, не против. У меня получилось реализовать только "автовыбор по убыванию"(не знаю как назвать), он будет только уменьшать шаг и если я введу уже изначально маленький шаг ,то автовыбор его не изменит и программа будет считать с ним.
0
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
0 / 0 / 0
Регистрация: 30.10.2016
Сообщений: 8
12.12.2016, 14:11  [ТС] 7
Извините , за некоторую назойливость , но не могли бы вы проверить мою программу ? Я прикрутил туда довольно топорный автовыбор шага ( на основе текста учебника ,еще я прочитал статью там использовался try_stepper, который я не очень понимаю , как реализовать , но его смысл понятен). Графики в особых точках соответствуют теории , но с погрешностями все равно беда.

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


C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
#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);
}
0
1481 / 1198 / 819
Регистрация: 29.02.2016
Сообщений: 3,579
12.12.2016, 15:41 8
Luckyscissors, выше я вам дал решение примерно той же задачи но другим методом, сравните с ней, подставте оттуда числа в свой код и сравните результаты
1
0 / 0 / 0
Регистрация: 30.10.2016
Сообщений: 8
12.12.2016, 23:15  [ТС] 9
На заданном в вашем решении промежутке времени различий между моим методом и вашим нету. Но как я понимаю , предоставленное решение методом Эйлера , сделано с целью увидеть аттрактор Лоренца ? В моей программе при r=28(когда
должен наблюдаться аттрактор) выходит какая-то мишура из точек (строю в Origin XYZ график для фазового портрета) , причем x(t), y(t) притягиваются к сложному множеству аттрактору.

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

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


C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
 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);
}
0
1481 / 1198 / 819
Регистрация: 29.02.2016
Сообщений: 3,579
13.12.2016, 07:55 10
можно продолжить сравнения, подставить теперь в метод Эйлера ваши данные и посмотреть что получится
1
0 / 0 / 0
Регистрация: 30.10.2016
Сообщений: 8
17.12.2016, 17:54  [ТС] 11
На больших промежутках времени даже предоставленное решение методом Эйлера теряет форму аттрактора Лоренца ( который по сути должен быть усточивым и принимать определенную форму при t->беск.). Может ли это быть ограничением метода ?
0
1481 / 1198 / 819
Регистрация: 29.02.2016
Сообщений: 3,579
17.12.2016, 18:03 12
метод Эйлера первого порядка, может быть такое, мне кажется
1
0 / 0 / 0
Регистрация: 30.10.2016
Сообщений: 8
18.12.2016, 23:07  [ТС] 13
Извините , тут обнаружилась еще одна ошибка. В готовом решении и моем , при r=19 наблюдаются неустойчивые фокусы(аттрактор Лоренца) , хотя при данном r должны получаться устойчивые фокусы.Аттрактор появляется при пороговом r>24,7. Исключая ошибки в теории , остаются только ошибки в коде. Вы можете предположить , что может быть не так ?
0
1481 / 1198 / 819
Регистрация: 29.02.2016
Сообщений: 3,579
19.12.2016, 09:14 14
попробуйте проверить это в еще какой то реализации метода, например в этой
https://www.codeproject.com/Ar... -equations
1
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
19.12.2016, 09:14
Помогаю со студенческими работами здесь

Метод рунге-кутта 3 порядка
дана функция d(y(x))/dx=e^x-2y(x) Нач. условия y(0)=e Код#include&lt;stdio.h&gt; #include&lt;conio.h&gt;...

Метод Рунге-Кутта четвертого порядка
Доброго времени суток. Задание по выч. мату. координаты х рассчитываются просто по шагу, но у...

Метод Рунге-Кутта второго порядка
В общем есть задание. Задание к-е на скрине: Ток в электрической цепи описывается...

Метод Рунге-Кутта 4 порядка (исправить код)
Всем привет. Нужно решить пример методом Рунге-Кутта 4-го порядка точности. Пример: y=1/2*x*y...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
14
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin
Copyright ©2000 - 2022, CyberForum.ru