Форум программистов, компьютерный форум, киберфорум
Наши страницы
Численные методы
Войти
Регистрация
Восстановить пароль
 
SVD102
0 / 0 / 1
Регистрация: 12.10.2015
Сообщений: 207
1

Не понятен метод RKF45

22.04.2016, 08:07. Просмотров 286. Ответов 7
Метки нет (Все метки)

Пишу программу по решению системы 2 уравнений 1 порядка методом RKF45. Не могу понять, как рассчитывается шаг на каждой итерации, и сам метод я не до конца понял. Мало информации в интернете, да и если есть то на англ. языке.
0
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
22.04.2016, 08:07
Ответы с готовыми решениями:

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

Не понятен метод связаный с boolean
мне не понятен полностью метод if(b) System.out.println("Эта инструкция...

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

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

не понятен алгоритм
Вобщем данный сценарий отбражает сообщения, есть два вида сообщений...

7
mathidiot
Эксперт по математике/физике
2888 / 2518 / 1105
Регистрация: 14.01.2014
Сообщений: 5,416
22.04.2016, 08:36 2
Это называется метод Рунге-Кутта-Фельберга.
Оптимальный шаг дается формулой (31)
0
Вложения
Тип файла: pdf RungeKuttaFehlbergProof.pdf (157.4 Кб, 3 просмотров)
SVD102
0 / 0 / 1
Регистрация: 12.10.2015
Сообщений: 207
22.04.2016, 08:52  [ТС] 3
Я уже нашел эту информацию, и не понял чему равен tol и чему равен h0 скажем, который мы задаем в начале, для первого шага
0
mathidiot
Эксперт по математике/физике
2888 / 2518 / 1105
Регистрация: 14.01.2014
Сообщений: 5,416
22.04.2016, 09:03 4
tol - это машинная точность, для первого шага выбирается любое разумное маленькое число, например: 0,001. В известных программах (Mathcad, Matlab) это значения по умолчанию
0
SVD102
0 / 0 / 1
Регистрация: 12.10.2015
Сообщений: 207
22.04.2016, 09:07  [ТС] 5
А чему равен h0. И еще про сам метод, то есть я посчитал новые значения Y и Z , далее что делать, можете алгоритм написать пожалуйста.
0
mathidiot
Эксперт по математике/физике
2888 / 2518 / 1105
Регистрация: 14.01.2014
Сообщений: 5,416
22.04.2016, 09:14 6
Прочитайте статью от начала до конца - там все подробно изложено! Повторять её содержимое здесь нет никакого смысла!
0
SVD102
0 / 0 / 1
Регистрация: 12.10.2015
Сообщений: 207
22.04.2016, 13:07  [ТС] 7
То есть получается посчитали Y и Z, сравнили их, если не удовлетворяют погрешности заданной то шаг уменьшается, а если удовлетворяют то шаг увеличивается. А какой первоначальный шаг брать, и это имеет значение какой шаг мы выберем?

Добавлено через 6 минут
Получается на каждой итерации h[i+1] = h[i]*s; Так?

Добавлено через 3 часа 20 минут
Вот код, но где - то здесь ошибка.

Кликните здесь для просмотра всего текста
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(fabs(y1 - x1) == eps/10.0||fabs(y1 - x1) < eps/10.0)
        {
            //h++;
            h = h*s(h,x1,y1);
            t = t + h;
 
        }
        else 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
        {
            //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
Том Ардер
Модератор
Эксперт по математике/физике
3835 / 2447 / 327
Регистрация: 15.06.2009
Сообщений: 4,473
22.04.2016, 13:52 8
SVD102, этот раздел - для обсуждения математических вопросов. Код и ошибки в нём следует обсуждать в соответствующем месте, напр. Метод RKF45?
0
22.04.2016, 13:52
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
22.04.2016, 13:52

не понятен код
Ребят не понятен. как бд подключать и т.п. unit Unit1; {главная...

Не понятен синтаксис...
WNDCLASS wc;//Объявление объекта оконного класса typedef struct...

Не понятен код
На практике приходится изучать COM и на одном из сайтов процессе обучения...


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

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

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Рейтинг@Mail.ru