Форум программистов, компьютерный форум, киберфорум
С++ для начинающих
Форум программистов и сисадминов Киберфорум > Форум программистов > C/C++ > C++ > С++ для начинающих
Войти
Регистрация
Восстановить пароль
Карта форума Темы раздела Блоги Сообщество Поиск Заказать работу  
 
Рейтинг 4.62/21: Рейтинг темы: голосов - 21, средняя оценка - 4.62
1 / 1 / 1
Регистрация: 12.10.2015
Сообщений: 207
1

Метод RKF45?

22.04.2016, 08:10. Показов 3981. Ответов 6
Метки нет (Все метки)
Author24 — интернет-сервис помощи студентам
Пишу программу по решению системы 2 уравнений 1 порядка методом RKF45. Не могу понять, как рассчитывается шаг на каждой итерации, и сам метод я не до конца понял. Мало информации в интернете, да и если есть то на англ. языке. Может подскажет кто-нибудь.
0
Лучшие ответы (1)
Сообщение: #9056385
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
 22.04.2016, 08:10
Ответы с готовыми решениями:

СЛАУ. Метод обратной матрицы, метод Гаусса, метод Крамера, метод Зейделя
Помогите ребят. Не могу построить алгоритмы для этих методов Язык C++

Метод медиан из трех элементов VS улучшенный быстрый метод сортировки(метод Бентли-Макилроя)
Здравствуйте! Дали весьма интересное задание. Сравнить два вышеуказанных метода сортировки для...

Мой код - метод бисекции, метод секущих (метод хорд)
Всем привет!!! Изучаем в институте С++. Сделал код, и там, и там одна и та же проблема - при любых...

Не понятен метод RKF45
Пишу программу по решению системы 2 уравнений 1 порядка методом RKF45. Не могу понять, как...

6
1550 / 875 / 179
Регистрация: 05.12.2015
Сообщений: 2,555
22.04.2016, 11:10 2
Цитата Сообщение от SVD102 Посмотреть сообщение
сам метод я не до конца понял
Для какого-то шага считаем рунге-кутту 4-го и 5-го порядков. Если между ними разница больше, чем требуется точность - уменьшаем шаг.
0
1 / 1 / 1
Регистрация: 12.10.2015
Сообщений: 207
22.04.2016, 12:00  [ТС] 3
И по новой пересчитываем эту итерацию с другим шагом уже ?
0
1550 / 875 / 179
Регистрация: 05.12.2015
Сообщений: 2,555
22.04.2016, 12:46 4
Цитата Сообщение от SVD102 Посмотреть сообщение
И по новой пересчитываем эту итерацию с другим шагом уже ?
Да, если в точность не влазит. В описании алгоритма должна быть формула вычисления шага. (Его же потом надо обратно увеличить)
0
1 / 1 / 1
Регистрация: 12.10.2015
Сообщений: 207
22.04.2016, 13:02  [ТС] 5
Вот мой код, но что -то не то выдает. Должно быть так: x1 =-0.90; x2 =0.47 (при t =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
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

#include <stdio.h>
#include <iostream>
#include <math.h>
using namespace std;
 
double const tol = 0.001;
 
double f1(double x1, double x2, double t)// 
{
    return(-52*x1-100*x2+exp(-t));               
}
 
double f2(double x1, double t)//
{
    return(x1+sin(t));               
} 
double s(double h, double x, double y)//
{
   return pow((tol*h)/(2*(fabs(y-x))),1.0/4.0);
}
 
int main()
{
    setlocale(LC_ALL,"rus");
   
    double t, s1,s2, h,tol, eps, E1, E2, k1, k2, k3, k4, k5, k6, y1, y2, m1, m2, m3, m4, m5, m6, x1, x2, x11,x22,y11,y22;
    //cout << "  "<<endl;
    //cout<<"h = ";//
    //cin >> h;
    
    t = 0; // 
    x1 = 1; // 
 
    x2 = 0; // 
    h = 0.001;   //2 or 3
    eps = 0.0001;//5 or 4
    int i = 0;
    y1 = x1;
    y2 = x2; 
    //s1 = s(h,x1,y1);
    //s2 = s(h,x2,y2);
    //while(E1 > eps || E2 > eps)
    while(true)
    {
 
        k1 = h*f1(x1, x2, t);
        m1 = h*f2(x1, t);
        
        k2 = h*f1(x1 + k1/4.0, x2 + m1/4.0, t + h/4.0);
        m2 = h*f2(x1 + k1/4.0, t + h/4.0);
 
        k3 = h*f1(x1 + (3.0*k1)/32.0 + (9.0*k2)/32.0, x2 + (3.0*m1)/32.0 + (9.0*m2)/32.0, t + (3.0*h)/8.0);
        m3 = h*f2(x1 + (3.0*k1)/32.0 + (9.0*k2)/32.0, t + (3.0*h)/8.0);
 
        k4 = h*f1(x1 + (1932.0*k1)/2197.0 - (7200.0*k2)/2197.0 + (7296.0*k3)/2197.0,
                  x2 + (1932.0*m1)/2197.0 - (7200.0*m2)/2197.0 + (7296.0*m3)/2197.0, t + (12.0*h)/13.0);
        m4 = h*f2(x1 + (1932.0*k1)/2197.0 - (7200.0*k2)/2197.0 + (7296.0*k3)/2197.0, t + (12.0*h)/13.0);
        
        k5 = h*f1(x1 + (439.0*k1)/216.0 - 8.0*k2 + (3680.0*k3)/513.0 - (815*k4)/4104.0,
                  x2 + (439.0*m1)/216.0 - 8.0*m2 + (3680.0*m3)/513.0 - (815*m4)/4104.0, t + h);
        m5 = h*f2(x1 + (439.0*k1)/216.0 - 8.0*k2 + (3680.0*k3)/513.0 - (815*k4)/4104.0, t + h);
 
        k6 = h*f1(x1 - (8.0*k1)/27.0 + 2*k2 - (3544.0*k3)/2565.0 + (1895.0*k4)/4104 - (11.0*k5)/40.0,
                  x2 - (8.0*m1)/27.0 + 2*m2 - (3544.0*m3)/2565.0 + (1895.0*m4)/4104 - (11.0*m5)/40.0,
                  t + h/2.0);
        m6 = h*f2(x1 - (8.0*k1)/27.0 + 2*k2 - (3544.0*k3)/2565.0 + (1895.0*k4)/4104 - (11.0*k5)/40.0, t + h/2.0);
 
        x11 = x1;
        x22 = x2;
        
        //до 5 
        x1 = x1 + (25.0*k1)/216.0 + (1408.0*k3)/2565.0 + (2197.0*k4)/4101.0 - k5/7.0;
        x2 = x2 + (25.0*m1)/216.0 + (1408.0*m3)/2565.0 + (2197.0*m4)/4101.0 - m5/7.0; 
        
        //cout<<"t = "<< t << " x1 = "<< x1 << " x2 = " << x2 << endl;
        //до 6
        y11 = x11;
        y22 = x22;
 
        y1 = y1 + (16.0*k1)/135.0 + (6656.0*k3)/12825.0 + (28561*k4)/56430.0 - (9*k5)/50.0 + (2.0*k6)/55.0; //???
        y2 = y2 + (16.0*m1)/135.0 + (6656.0*m3)/12825.0 + (28561*m4)/56430.0 - (9*m5)/50.0 + (2.0*m6)/55.0; //???
        //if(fabs(x1-x11)< 0.0001||fabs(x2-x22)< 0.0001)
        //{
        //  break;
        //}
        if((eps<=fabs(y1 - x1) && fabs(y1 - x1) < eps*10.0) && (eps<=fabs(y1 - x1)&&fabs(y1 - x1) < eps*10.0))
        {
            cout<<"Ответ:"<<endl;
            cout<<"t = "<< t << " x1 = "<< x1<< " y1 = "<< y1<< " x2 = "<< x2<< " y2 = "<< y2<<endl;
            break;
        }
        else if(fabs(y1 - x1) == eps/10.0||fabs(y1 - x1) < eps/10.0)
        {
            //h++;
            h = h*s(h,x1,y1);
            t = t + h;
 
        }
        else
        {
            //h--;
            h = h*s(h,x1,y1);
            //t = t + h;
            x1 = x11;
            x2 = x22;
            y1 = y11;
            y2 = y22;
        }
        //cмотрим шаг, если погрешность не устраивает . то меняем шаг и тогда удаляем  последнюю итерацию(вроде)
        //если погрешность устраивает то шаг оставляем и продолжаем.
        //E1 = ();
        //cout<<"E1 = "<<E1<<"  E2 = "<<E2;
        
    }
    
    //cout<<"t = "<< t << " x1 = "<< x1 << " x2 = " << x2 << endl;
    
    system("pause");
    return 0;
}
0
1550 / 875 / 179
Регистрация: 05.12.2015
Сообщений: 2,555
22.04.2016, 15:52 6
Лучший ответ Сообщение было отмечено SVD102 как решение

Решение

SVD102, Исправил. Увеличение шага, пределы и т.д. сами прикрутите.
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

#include <stdio.h>
#include <iostream>
#include <math.h>
using namespace std;
 
double const tol = 0.0001;
 
double f1(double x1, double x2, double t)// 
{
    return(-52 * x1 - 100 * x2 + exp(-t));
}
 
double f2(double x1, double t)//
{
    return(x1 + sin(t));
}
double s(double h, double x, double y)//
{
    return pow((tol*h) / (2 * (fabs(y - x))), 1.0 / 4.0);
}
 
int main()
{
    setlocale(LC_ALL, "rus");
 
    double t, s1, s2, h, tol, eps, E1, E2, k1, k2, k3, k4, k5, k6, y1, y2, m1, m2, m3, m4, m5, m6, x1, x2, x11, x22, y11, y22;
    //cout << "  "<<endl;
    //cout<<"h = ";//
    //cin >> h;
 
    t = 0; // 
    x1 = 1; // 
 
    x2 = 0; // 
    h = 0.001;   //2 or 3
    eps = 0.0001;//5 or 4
    int i = 0;
    y1 = x1;
    y2 = x2;
    //s1 = s(h,x1,y1);
    //s2 = s(h,x2,y2);
    //while(E1 > eps || E2 > eps)
    while (true)
    {
 
        k1 = h*f1(x1, x2, t);
        m1 = h*f2(x1, t);
 
        k2 = h*f1(x1 + k1 / 4.0, x2 + m1 / 4.0, t + h / 4.0);
        m2 = h*f2(x1 + k1 / 4.0, t + h / 4.0);
 
        k3 = h*f1(x1 + (3.0*k1) / 32.0 + (9.0*k2) / 32.0, x2 + (3.0*m1) / 32.0 + (9.0*m2) / 32.0, t + (3.0*h) / 8.0);
        m3 = h*f2(x1 + (3.0*k1) / 32.0 + (9.0*k2) / 32.0, t + (3.0*h) / 8.0);
 
        k4 = h*f1(x1 + (1932.0*k1) / 2197.0 - (7200.0*k2) / 2197.0 + (7296.0*k3) / 2197.0,
            x2 + (1932.0*m1) / 2197.0 - (7200.0*m2) / 2197.0 + (7296.0*m3) / 2197.0, t + (12.0*h) / 13.0);
        m4 = h*f2(x1 + (1932.0*k1) / 2197.0 - (7200.0*k2) / 2197.0 + (7296.0*k3) / 2197.0, t + (12.0*h) / 13.0);
 
        k5 = h*f1(x1 + (439.0*k1) / 216.0 - 8.0*k2 + (3680.0*k3) / 513.0 - (815 * k4) / 4104.0,
            x2 + (439.0*m1) / 216.0 - 8.0*m2 + (3680.0*m3) / 513.0 - (815 * m4) / 4104.0, t + h);
        m5 = h*f2(x1 + (439.0*k1) / 216.0 - 8.0*k2 + (3680.0*k3) / 513.0 - (815 * k4) / 4104.0, t + h);
 
        k6 = h*f1(x1 - (8.0*k1) / 27.0 + 2 * k2 - (3544.0*k3) / 2565.0 + (1895.0*k4) / 4104 - (11.0*k5) / 40.0,
            x2 - (8.0*m1) / 27.0 + 2 * m2 - (3544.0*m3) / 2565.0 + (1895.0*m4) / 4104 - (11.0*m5) / 40.0,
            t + h / 2.0);
        m6 = h*f2(x1 - (8.0*k1) / 27.0 + 2 * k2 - (3544.0*k3) / 2565.0 + (1895.0*k4) / 4104 - (11.0*k5) / 40.0, t + h / 2.0);
 
        //до 5 
        double delta51 = (25.0*k1) / 216.0 + (1408.0*k3) / 2565.0 + (2197.0*k4) / 4101.0 - k5 / 7.0;
        double delta52 = (25.0*m1) / 216.0 + (1408.0*m3) / 2565.0 + (2197.0*m4) / 4101.0 - m5 / 7.0;
 
        double delta61 = (16.0*k1) / 135.0 + (6656.0*k3) / 12825.0 + (28561 * k4) / 56430.0 - (9 * k5) / 50.0 + (2.0*k6) / 55.0; //???
        double delta62 = (16.0*m1) / 135.0 + (6656.0*m3) / 12825.0 + (28561 * m4) / 56430.0 - (9 * m5) / 50.0 + (2.0*m6) / 55.0; //???
 
        if (fabs(delta61 - delta51) > eps || fabs(delta62 - delta52) > eps) {
            h *= s(h, delta51, delta61);
            continue;
        }
        x1 += delta61;
        x2 += delta62;
        t += h;
        
        if (2 - t < h)
            h = 2 - t;
        if (t >= 2)
            break;
 
        //cмотрим шаг, если погрешность не устраивает . то меняем шаг и тогда удаляем  последнюю итерацию(вроде)
        //если погрешность устраивает то шаг оставляем и продолжаем.
        //E1 = ();
        //cout<<"E1 = "<<E1<<"  E2 = "<<E2;
 
    }
 
    cout<<"t = "<< t << " x1 = "<< x1 << " x2 = " << x2 << endl;
 
    system("pause");
    return 0;
}
Опять вы какие-то x11, y22 напихали.

На будущее: китайский код - он такой китайский.
Перепишите это все через произведение матрицы на вектор.
0
1 / 1 / 1
Регистрация: 12.10.2015
Сообщений: 207
25.04.2016, 08:03  [ТС] 7
C++
1
2

if (2 - t < h)
h = 2 - t;
Откуда это условие, зачем оно нужно?
0
25.04.2016, 08:03
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
 25.04.2016, 08:03
Помогаю со студенческими работами здесь

Runge-Kutta-Fehlberg (RKF45)
Здравствуйте, мне нужно решить численным методом Runge-Kutta-Fehlberg (RKF45) дифференциальное...

Не пойму как исправить ошыбку (RKF45)
Код решает дифференциальное уравнение численным метдом &quot;RKF45&quot;, но не пойму в чем ошыбка. Буду рад...

Написать три алгоритма решения СЛАУ: Метод прогонки, метод квадратных корней, метод вращений
Начал писать курсовую. Нужно написать три алгоритма решения СЛАУ: прогонки, квадратных корней,...

Метод хорд, метод касательных, метод половинного деления
Ребят помогите пожалуйста. Задали тему курсовой &quot;Решение нелинейных уравнений (метод хорд, метод...


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

Или воспользуйтесь поиском по форуму:
7
Ответ Создать тему
КиберФорум - форум программистов, компьютерный форум, программирование
Поддержать - Реклама - Обратная связь
Powered by vBulletin
Copyright ©2000 - 2024, CyberForum.ru