метод рунге-кутты-мерсона

@mrniws · Регистрация: 16.01.2010

Студворк — интернет-сервис помощи студентам

нужно написать программу для решения систем из n дифф уравнений методом рунге-кутты-мерсона(он же пятиэтапный метод рунге кутта 4 порядка).

вот кусок кода,который считает для одного уравнения.написан не идеально,только учусь..вот не могу понять как сделать для системы..уже сижу дня три и уже голову сломал..есть еще куски,пробовал написать для двух уравнений,но такое не принимают...
спасите меня...хоть какой-нибудь идейкой...просто исчерпал все свои знания.
у нас есть система вида
y1'=a1*y1+b1*y2+..+N1*yn
y2'=a2*y1+b2*y2+..+N2*yn
.....................................
yn'=an*y1+bn*y2+..+Nn*yn
и есть начальные условия,как для задачи коши

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
do
    {
        k1=f(x[i],y[i]);
        k2=f(x[i]+h/3,y[i]+h/3*k1);
        k3=f(x[i]+h/3,y[i]+h/6*k1+h/6*k2);
        k4=f(x[i]+h/2,y[i]+h/8*k1+3*h/8*k2);
        k5=f(x[i]+h,y[i]+h/2*k1-3*h/2*k3+2*h*k4);
        yx=y[i]+h/2*(k1-3*k3+4*k4);
        y=y[i]+h/6*(k1+4*k4+k5);
        r=0.2*abs(y-yx);
        if(r>e) h=h/2;
        else 
        {
            if((r<=e/64)&&(2*h<h0)) h=2*h;
            else h=h;
        }
        while(r>e)
        {
            k1=f(x[i],y[i]);
            k2=f(x[i]+h/3,y[i]+h/3*k1);
            k3=f(x[i]+h/3,y[i]+h/6*k1+h/6*k2);
            k4=f(xi]+h/2,y[i]+h/8*k1+3*h/8*k2);
            k5=f(x[i]+h,y[i]+h/2*k1-3*h/2*k3+2*h*k4);
            yx=y[i]+h/2*(k1-3*k3+4*k4);
            y=y[i]+h/6*(k1+4*k4+k5);
            r=0.2*abs(y-yx);
            if(r>e) h=h/2;
            else 
            {
                if((r<=e/64)&&(2*h<h0)) h=2*h;
                else h=h;
            }
        }
 
        y[i+1]=y;
        i=i+1;
        x[i]=x[i-1]+h;
        
    }
    while (x[i-1]<=xn);

@Hagacure · 02.01.2011, 02:57

Вот вольный перевод на Си из паскаля паскалевская версия взята из книги Мудрова "Численные методы ПЭКМ"

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
double X,X9,H, EH,EPS,H0,P,D,X1,Y[8];
//------------------------------------------------------ ïðîèçâîäíûå *>
void RP(double X,double *Y,double *F)
{
   F[0] = Y[1];
   F[1] = P * (1.0 - Y[0]*Y[0])*Y[1]-Y[0];
}
//------------------------------------------------ ÌÅÒÎÄ ÐÓÍÃÅ-ÊÓÒÒÛ-ÌÅÐÑÎÍÀ *>
void RKM(int N,double *X,double *H,double *E, double *Y)
{
        int I;
        double H3,H4,R = 0.0,A,Z[8],F0[8],F[8],K1[8],K3[8];
        RP(*X,&Y[0],&F[0]);
        for (I=0;I<N ;I++)
           Z[I] = Y[I];
        do
        {
                H3 = *H / 3.0;
                H4 = 4 * H3;
                for (I=0; I<N; I++)
                {
                        K1[I] = H3 * F0[I];
                        Y[I] = Z[I] + K1[I];
                }
                RP(*X + H3,&Y[0],&F[0]);
                for (I=0; I<N; I++)
                   Y[I] = Z[I] + (K1[I] + H3*F[I])/2;
                RP(*X + H3,&Y[0],&F[0]);
                for (I=0; I<N; I++)
                {
                        K3[I] = *H * F[I];
                        Y[I] = Z[I] + 0.375 * (K1[I] + K3[I]);
                }
                RP((*X + *H/2),&Y[0],&F[0]);
                for (I=0; I<N; I++)
                {
                        K1[I] += H4 * F[I];
                        Y[I] = Z[I] + 1.5 * (K1[I] + K3[I]);
                }
                RP(*X + *H/2,&Y[0],&F[0]);
                for (I=0; I<N; I++)
                {
                        K1[I] += H4 * F[I];
                        Y[I] = Z[I] + 1.5 * (K1[I] - K3[I]);
                }
                RP(*X + *H,&Y[0],&F[0]);
                for (I=0; I<N; I++)
                {
                        A = H3 * F[I];
                        Y[I] = Z[I] + (K1[I] + A)/2;
                        A = 2 * K1[I] - 3.0 * K3[I] - A;
                        if ( Y[I] != 0.0)
                                A /= Y[I];
                        if ( abs(A) > R )
                                R=abs(A);
                }
                *H = *H/2;
        }
        while(R<*E);
        *H *= 2;
        *X += *H;
        if (32*R < *E)
                *H *= 2;
}
void __fastcall TForm1::run() 
{
//входящие данные
        X=Edit1->Text.ToDouble();
        X9=Edit2->Text.ToDouble();
        H=Edit3->Text.ToDouble();
        EH=Edit4->Text.ToDouble();
        EPS=Edit5->Text.ToDouble();
        H0=Edit6->Text.ToDouble();
        Y[0]=Edit7->Text.ToDouble();
        Y[1]=Edit8->Text.ToDouble();
        P=Edit9->Text.ToDouble();
//входящие данные
        X1 = X;
        RKM(2,&X,&H0,&EPS,Y);
        H0 =  Sign (H);
        D  = X1 + H - X;
        if (D < 0.0 == (H>0.0))
                H0 = D;
}

@mrniws · 06.01.2011, 17:52 **[ТС]**

спасибо=)

@Boris) · 07.05.2011, 18:18

А есть способ решения (делфи или паскаль или си) для 1-го диф. уравнения n-го порядка y^n=f(y^(n-1),...,t) при заданых начальных условиях ?

Новые блоги и статьи Все статьи Все блоги /
PhpStorm 2025.3: WSL Terminal всегда стартует в ~ and_y87 14.12.2025 PhpStorm 2025. 3: WSL Terminal всегда стартует в ~ (home), игнорируя директорию проекта Симптом: После обновления до PhpStorm 2025. 3 встроенный терминал WSL открывается в домашней директории. . .	Access VikBal 11.12.2025 Помогите пожалуйста !! Как объединить 2 одинаковые БД Access с разными данными.	Новый ноутбук volvo 07.12.2025 Всем привет. По скидке в "черную пятницу" взял себе новый ноутбук Lenovo ThinkBook 16 G7 на Амазоне: Ryzen 5 7533HS 64 Gb DDR5 1Tb NVMe 16" Full HD Display Win11 Pro	Музыка, написанная Искусственным Интеллектом volvo 04.12.2025 Всем привет. Некоторое время назад меня заинтересовало, что уже умеет ИИ в плане написания музыки для песен, и, собственно, исполнения этих самых песен. Стихов у нас много, уже вышли 4 книги, еще 3. . .	От async/await к виртуальным потокам в Python IndentationError 23.11.2025 Армин Ронахер поставил под сомнение async/ await. Создатель Flask заявляет: цветные функции - провал, виртуальные потоки - решение. Не threading-динозавры, а новое поколение лёгких потоков. Откат?. . .
Поиск "дружественных имён" СОМ портов Argus19 22.11.2025 Поиск "дружественных имён" СОМ портов На странице: https:/ / norseev. ru/ 2018/ 01/ 04/ comportlist_windows/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .	Сколько Государство потратило денег на меня, обеспечивая инсулином. Programma_Boinc 20.11.2025 Сколько Государство потратило денег на меня, обеспечивая инсулином. Вот решила сделать интересный приблизительный подсчет, сколько государство потратило на меня денег на покупку инсулинов. . . .	Ломающие изменения в C#.NStar Alpha Etyuhibosecyu 20.11.2025 Уже можно не только тестировать, но и пользоваться C#. NStar - писать оконные приложения, содержащие надписи, кнопки, текстовые поля и даже изображения, например, моя игра "Три в ряд" написана на этом. . .	Мысли в слух kumehtar 18.11.2025 Кстати, совсем недавно имел разговор на тему медитаций с людьми. И обнаружил, что они вообще не понимают что такое медитация и зачем она нужна. Самые базовые вещи. Для них это - когда просто люди. . .	Создание Single Page Application на фреймах krapotkin 16.11.2025 Статья исключительно для начинающих. Подходы оригинальностью не блещут. В век Веб все очень привыкли к дизайну Single-Page-Application . Быстренько разберем подход "на фреймах". Мы делаем одну. . .

@mrniws 0 / 0 / 0 Регистрация: 16.01.2010 Сообщений: 3
	06.01.2011, 17:52 [ТС]
	спасибо=) 0

@Boris) 0 / 0 / 0 Регистрация: 06.05.2011 Сообщений: 12
	07.05.2011, 18:18
	А есть способ решения (делфи или паскаль или си) для 1-го диф. уравнения n-го порядка y^n=f(y^(n-1),...,t) при заданых начальных условиях ? 0