Какой способ лучше для решения данной системы дифуров?

@-Sky- · Регистрация: 22.09.2010

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

Привет!)
У меня проблемка, не могу решить системы дифуров численно.

Система:
dx/dt=a1*x-b1*x*y
dy/dt=-a2*x+b2*x*y
где a1,a2,b1,b2 произвольные константы (пусть они равны еденице), начальные параметры и отрезок дифференцирования тоже не суть важно

Дело даже не в том что я не могу написать код, а в том что происходит катастрофическое накопление погрешностей, если я использую например мотод Ранге-Кутта 4ого порядка или метод прогноз-коррекция. Возможно дело в том что общее решение включает в себя экспоненты, ... хотя какой дифур без них обходится...
С методом Милна (на систему я его так и не обобщил) у меня вообще что то странное, происходит накопление погрешностей даже на ровном месте.... может я не правильно код написал???

Может кто укажет на ошибку но вроде все правильно (в функциях f1 и f2 задана другая система, я ее использовал для проверки правельности кода):
Метод Рунге-Кутта 4ого порядка:

Delphi
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
program S2R4;
 
{$APPTYPE CONSOLE}
 
uses
  SysUtils;
 
type tA=Array of Extended;
 
function f1(t,x,y:Extended):Extended;
var a1,b1:Extended;
begin
  a1:=1;
  b1:=1;
  result:=y;
end;
 
function f2(t,x,y:Extended):Extended;
var a2,b2:Extended;
begin
  a2:=1;
  b2:=1;
  result:=(x/t-y)/t-x;
end;
 
var h,a,b,dx,dy:Extended;
    k1,k2,k3,k4,l1,l2,l3,l4:Extended;
    N,i:Integer;
    T,X,Y:tA;
    fT,fX,fY:TextFile;
begin
 
  assignFile(fT,'outputT.txt');
  rewrite(fT);
  assignFile(fX,'outputX.txt');
  rewrite(fX);
  assignFile(fY,'outputY.txt');
  rewrite(fY);
 
  a:=1;
  b:=10;
  h:=0.01;
 
  N:=trunc(abs(b-a)/h);
 
  setLength(T,N+1);
  setLength(X,N+1);
  setLength(Y,N+1);
 
  T[0]:=a;
  X[0]:=0.1;
  Y[0]:=0.5;
 
  for i:=0 to N do
  begin
 
    writeln(fT,T[i]);
    writeln(fX,X[i]);
    writeln(fY,Y[i]);
 
    k1:=h*f1(T[i],X[i],Y[i]);
    l1:=h*f2(T[i],X[i],Y[i]);
 
    k2:=h*f1(T[i]+h/2,X[i]+k1/2,Y[i]+l1/2);
    l2:=h*f2(T[i]+h/2,X[i]+k1/2,Y[i]+l1/2);
 
    k3:=h*f1(T[i]+h/2,X[i]+k2/2,Y[i]+l2/2);
    l3:=h*f2(T[i]+h/2,X[i]+k2/2,Y[i]+l2/2);
 
    k4:=h*f1(T[i]+h,X[i]+k3,Y[i]+l3);
    l4:=h*f2(T[i]+h,X[i]+k3,Y[i]+l3);
 
    dx:=(k1+2*k2+2*k3+k4)/6;
    dy:=(l1+2*l2+2*l3+l4)/6;
 
    T[i+1]:=T[i]+h;
    X[i+1]:=X[i]+dx;
    Y[i+1]:=Y[i]+dy;
 
  end;
 
  closeFile(fT);
  closeFile(fX);
  closeFile(fY);
end.

Метод Прогноз-Коррекция:

Delphi
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
program PK;
 
{$APPTYPE CONSOLE}
 
uses
  SysUtils;
 
type tA=Array of Extended;
 
function f1(x,y,z:Extended):Extended;
var a1,b1:Extended;
begin
  a1:=1;
  b1:=1;
  result:=z;
end;
 
function f2(x,y,z:Extended):Extended;
var a2,b2:Extended;
begin
  a2:=1;
  b2:=1;
  result:=(y/x-z)/x-y;
end;
 
var h,a,b,Yp12,Zp12:Extended;
    N,i:Integer;
    X,Y,Yp,Z,Zp:tA;
    fX,fY,fZ:TextFile;
begin
 
  assignFile(fX,'outputX.txt');
  rewrite(fX);
  assignFile(fY,'outputY.txt');
  rewrite(fY);
  assignFile(fZ,'outputZ.txt');
  rewrite(fZ);
 
  a:=1;
  b:=10;
  h:=0.01;
 
  N:=trunc(abs(b-a)/h);
 
  setLength(X,N+1);
  setLength(Y,N+1);
  setLength(Yp,N+1);
  setLength(Z,N+1);
  setLength(Zp,N+1);
 
  X[0]:=a;
  Y[0]:=0.1;
  Yp[0]:=0.1;
  Z[0]:=0.5;
  Zp[0]:=0.5;
 
  Yp12:=Yp[0]+h*f1(X[0],Yp[0],Zp[0])/2;
  Zp12:=Zp[0]+h*f2(X[0],Yp[0],Zp[0])/2;
 
  Yp[1]:=Yp[0]+h*f1(a+h/2,Yp12,Zp12);
  Zp[1]:=Zp[0]+h*f2(a+h/2,Yp12,Zp12);
 
  for i:=0 to N do
  begin
 
    writeln(fX,X[i]);
    writeln(fY,Y[i]);
    writeln(fZ,Z[i]);
 
    X[i+1]:=X[i]+h;
    Y[i+1]:=Y[i]+h*(f1(X[i],Yp[i],Zp[i])+f1(X[i+1],Yp[i+1],Zp[i+1]))/2;
    Yp[i+2]:=Y[i]+2*h*f1(X[i+1],Y[i+1],Z[i+1]);
 
    Z[i+1]:=Z[i]+h*(f2(X[i],Yp[i],Zp[i])+f2(X[i+1],Yp[i+1],Zp[i+1]))/2;
    Zp[i+2]:=Z[i]+2*h*f2(X[i+1],Y[i+1],Z[i+1]);
 
  end;
 
  closeFile(fX);
  closeFile(fY);
  closeFile(fZ);
 
end.

а вот с методом Милна совсем что то странное... может не правильная реализация? Я пока не стал писать систему, рассмотрел только одномерный случай.

Метод Милна:

Delphi
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
program S1Miln;
 
{$APPTYPE CONSOLE}
 
uses
  SysUtils;
 
type tA=Array of Extended;
 
function f1(t,x,y:Extended):Extended;
var a1,b1:Extended;
begin
  a1:=1;
  b1:=1;
  result:=-2*t*sqr(x);
end;
 
function f2(t,x,y:Extended):Extended;
var a2,b2:Extended;
begin
  a2:=1;
  b2:=1;
  result:=-a2*x+b2*x*y;
end;
 
var h,a,b,dx,dy:Extended;
    k1,k2,k3,k4,l1,l2,l3,l4:Extended;
    N,i:Integer;
    T,X,Y,Xp,Yp:tA;
    fT,fX,fY:TextFile;
begin
 
// на массивы Y и Yp не обращайте внимания т.к. тут решается ОДНОМЕРНАЯ СИСТЕМА
 
  assignFile(fT,'outputT.txt');
  rewrite(fT);
  assignFile(fX,'outputX.txt');
  rewrite(fX);
  assignFile(fY,'outputY.txt');
  rewrite(fY);
 
  a:=0;
  b:=1;
  h:=0.1;
 
  N:=trunc(abs(b-a)/h);
 
  setLength(T,N+1);
  setLength(X,N+1);
  setLength(Y,N+1);
  setLength(Xp,N+1);
  setLength(Yp,N+1);
 
  T[0]:=a;
  X[0]:=1;
  Y[0]:=0.5;
  Xp[0]:=1;
  Yp[0]:=0.5;
 
  for i:=0 to 3 do
  begin
 
    k1:=h*f1(T[i],Xp[i],Yp[i]);
    l1:=h*f2(T[i],Xp[i],Yp[i]);
 
    k2:=h*f1(T[i]+h/2,Xp[i]+k1/2,Yp[i]+l1/2);
    l2:=h*f2(T[i]+h/2,Xp[i]+k1/2,Yp[i]+l1/2);
 
    k3:=h*f1(T[i]+h/2,Xp[i]+k2/2,Yp[i]+l2/2);
    l3:=h*f2(T[i]+h/2,X[pi]+k2/2,Yp[i]+l2/2);
 
    k4:=h*f1(T[i]+h,Xp[i]+k3,Yp[i]+l3);
    l4:=h*f2(T[i]+h,Xp[i]+k3,Yp[i]+l3);
 
    dx:=(k1+2*k2+2*k3+k4)/6;
    dy:=(l1+2*l2+2*l3+l4)/6;
 
    T[i+1]:=T[i]+h;
    Xp[i+1]:=Xp[i]+dx;
    Yp[i+1]:=Yp[i]+dy;
 
  end;
 
//  Xp[1]:=0.9901;
//  Xp[2]:=0.9615;
//  Xp[3]:=0.9173;
 
 
  X[1]:=Xp[1];
 
  for i:=1 to 2 do
  begin
 
    X[i+1]:=X[i]+(h/3)*(f1(T[i-1],Xp[i-1],Yp[i-1])+4*f1(T[i],Xp[i],Yp[i])+f1(T[i+1],Xp[i+1],Yp[i+1]));
 
  end;
 
  for i:=3 to N do
  begin
 
    T[i+1]:=T[i]+h;
    X[i+1]:=X[i]+(h/3)*(f1(T[i-1],Xp[i-1],Yp[i-1])+4*f1(T[i],Xp[i],Yp[i])+f1(T[i+1],Xp[i+1],Yp[i+1]));
    Xp[i+1]:=X[i-3]+(4*h/3)*(2*f1(T[i-2],X[i-2],Y[i-2])-f1(T[i-1],X[i-1],Y[i-1])+2*f1(T[i],X[i],Y[i]));
 
  end;
 
  for i:=0 to N do
  begin
 
    writeln(fT,T[i]);
    writeln(fX,X[i]);
    writeln(fY,Y[i]);
 
  end;
 
 
  closeFile(fT);
  closeFile(fX);
  closeFile(fY);
end.

Помогите кто чем может

KuKu · 22.09.2010, 21:23

аналитически попробуй

, по размеру как раз выйдет как на один из методов, даже чуть меньше. ))

@-Sky- · 22.09.2010, 21:33 **[ТС]**

Так в том то все и дело, что аналитически не получится! там слишком кривая ерунда получается...

KuKu · 22.09.2010, 22:33

Шаг попробуй поуменьшать. f1 и f2 у тебя какие то мутные еще...

Добавлено через 13 минут
это про рунге-кутта говорю, про остальные не знаю

@-Sky- · 22.09.2010, 23:23 **[ТС]**

да ... там константы просто остались от старой системы, а так все нормально с ними

KuKu · 23.09.2010, 00:21

C++
1
result:=y;

вродь у рунге-кута решаем уравнение вида x'=f(x,t) и y=f(y,t) , что просто есть уравнения без дополнительных преобразований.

@-Sky- · 23.09.2010, 00:39 **[ТС]**

Я уже говорил, что в функциях f1 и f2 задана другая система, а именно:
x'=y
y'=(x/t-y)/t-x
я ее ввел для проверки на правельность кода
мне же нужно решить
x'=a1*x-b1*x*y
y'=-a2*x+b2*x*y
да и это можно записать как:
x'=f1(x,y)
y'=f2(x,y) от t моя система не зависит

Добавлено через 58 секунд
то есть система которую мне нужно решить

KuKu · 23.09.2010, 03:38

метод рунге-кутта решает уравнение вида x'=f1(x,y), когда x' - производная по y(т.е dx/dy), у тебя же это x'=f1(x,y), где x' производная по времени(от 3ей переменной).

@-Sky- · 23.09.2010, 04:04 **[ТС]**

ты ошибся, то что ты написал верно для одномерного случая, а у меня система из 2ух дифуров и у меня 3 переменных (пусть будут не t,x,y как у меня, а x,y,z), тогда в общем случае можно написать
dy/dx=f1(x,y,z)
dz/dx=f2(x,y,z)
тогда то и получается что ты написал

KuKu · 23.09.2010, 11:29

ну если так то извиняюсь

. Ошибка хоть как будет при многом количестве итерраций. Опять же можно уменьшать шаг или попробуй описать приближения более высоких порядков. Можно попробывать часть решить аналитически, связь между x и y находится легко аналитчески это точно. Дальше решать только одно уравнение, ошибку это должно уменьшить, но не думаю что существенно. Можешь поставить математические пакеты и пусть они попробуют решить аналитически ... не надо пренебрегать теорие в данном случае, в аналитическом нет ошибок вообще. На форуме есть раздел дифуров, там думаю математиков найдешь, которые тебе скажут точно есть ли решение у подобных уравнений.

@-Sky- · 24.09.2010, 20:35 **[ТС]**

Всем спасибо, я сам разобрался!

Новые блоги и статьи Все статьи Все блоги /
Инструменты COM: Сохранение данный из VARIANT в файл и загрузка из файла в VARIANT bedvit 28.01.2026 Сохранение базовых типов COM и массивов (одномерных или двухмерных) любой вложенности (деревья) в файл, с возможностью выбора алгоритмов сжатия и шифрования. Часть библиотеки BedvitCOM Использованы. . .	Загрузка PNG с альфа-каналом на SDL3 для Android: с помощью SDL_LoadPNG (без SDL3_image) 8Observer8 28.01.2026 Содержание блога SDL3 имеет собственные средства для загрузки и отображения PNG-файлов с альфа-каналом и базовой работы с ними. В этой инструкции используется функция SDL_LoadPNG(), которая. . .	Загрузка PNG с альфа-каналом на SDL3 для Android: с помощью SDL3_image 8Observer8 27.01.2026 Содержание блога SDL3_image - это библиотека для загрузки и работы с изображениями. Эта пошаговая инструкция покажет, как загрузить и вывести на экран смартфона картинку с альфа-каналом, то есть с. . .	Влияние грибов на сукцессию anaschu 26.01.2026 Бифуркационные изменения массы гриба происходят тогда, когда мы уменьшаем массу компоста в 10 раз, а скорость прироста биомассы уменьшаем в три раза. Скорость прироста биомассы может уменьшаться за. . .
Воспроизведение звукового файла с помощью SDL3_mixer при касании экрана Android 8Observer8 26.01.2026 Содержание блога SDL3_mixer - это библиотека я для воспроизведения аудио. В отличие от инструкции по добавлению текста код по проигрыванию звука уже содержится в шаблоне примера. Нужно только. . .	Установка Android SDK, NDK, JDK, CMake и т.д. 8Observer8 25.01.2026 Содержание блога Перейдите по ссылке: https:/ / developer. android. com/ studio и в самом низу страницы кликните по архиву "commandlinetools-win-xxxxxx_latest. zip" Извлеките архив и вы увидите. . .	Вывод текста со шрифтом TTF на Android с помощью библиотеки SDL3_ttf 8Observer8 25.01.2026 Содержание блога Если у вас не установлены Android SDK, NDK, JDK, и т. д. то сделайте это по следующей инструкции: Установка Android SDK, NDK, JDK, CMake и т. д. Сборка примера Скачайте. . .	Использование SDL3-callbacks вместо функции main() на Android, Desktop и WebAssembly 8Observer8 24.01.2026 Содержание блога Если вы откроете примеры для начинающих на официальном репозитории SDL3 в папке: examples, то вы увидите, что все примеры используют следующие четыре обязательные функции, а. . .

KuKu 1563 / 1041 / 94 Регистрация: 17.04.2009 Сообщений: 2,995
	22.09.2010, 21:23
	аналитически попробуй , по размеру как раз выйдет как на один из методов, даже чуть меньше. )) 0

@-Sky- 3 / 3 / 0 Регистрация: 22.09.2010 Сообщений: 18
	22.09.2010, 21:33 [ТС]
	Так в том то все и дело, что аналитически не получится! там слишком кривая ерунда получается... 0

KuKu 1563 / 1041 / 94 Регистрация: 17.04.2009 Сообщений: 2,995
	22.09.2010, 22:33
	Шаг попробуй поуменьшать. f1 и f2 у тебя какие то мутные еще... Добавлено через 13 минут это про рунге-кутта говорю, про остальные не знаю 0

@-Sky- 3 / 3 / 0 Регистрация: 22.09.2010 Сообщений: 18
	22.09.2010, 23:23 [ТС]
	да ... там константы просто остались от старой системы, а так все нормально с ними 0

@-Sky- 3 / 3 / 0 Регистрация: 22.09.2010 Сообщений: 18
	23.09.2010, 00:39 [ТС]
	Я уже говорил, что в функциях f1 и f2 задана другая система, а именно: x'=y y'=(x/t-y)/t-x я ее ввел для проверки на правельность кода мне же нужно решить x'=a1x-b1xy y'=-a2x+b2xy да и это можно записать как: x'=f1(x,y) y'=f2(x,y) от t моя система не зависит Добавлено через 58 секунд то есть система которую мне нужно решить 0

KuKu 1563 / 1041 / 94 Регистрация: 17.04.2009 Сообщений: 2,995
	23.09.2010, 03:38
	метод рунге-кутта решает уравнение вида x'=f1(x,y), когда x' - производная по y(т.е dx/dy), у тебя же это x'=f1(x,y), где x' производная по времени(от 3ей переменной). 0

@-Sky- 3 / 3 / 0 Регистрация: 22.09.2010 Сообщений: 18
	23.09.2010, 04:04 [ТС]
	ты ошибся, то что ты написал верно для одномерного случая, а у меня система из 2ух дифуров и у меня 3 переменных (пусть будут не t,x,y как у меня, а x,y,z), тогда в общем случае можно написать dy/dx=f1(x,y,z) dz/dx=f2(x,y,z) тогда то и получается что ты написал 0

KuKu 1563 / 1041 / 94 Регистрация: 17.04.2009 Сообщений: 2,995
	23.09.2010, 11:29
	ну если так то извиняюсь. Ошибка хоть как будет при многом количестве итерраций. Опять же можно уменьшать шаг или попробуй описать приближения более высоких порядков. Можно попробывать часть решить аналитически, связь между x и y находится легко аналитчески это точно. Дальше решать только одно уравнение, ошибку это должно уменьшить, но не думаю что существенно. Можешь поставить математические пакеты и пусть они попробуют решить аналитически ... не надо пренебрегать теорие в данном случае, в аналитическом нет ошибок вообще. На форуме есть раздел дифуров, там думаю математиков найдешь, которые тебе скажут точно есть ли решение у подобных уравнений. 0

@-Sky- 3 / 3 / 0 Регистрация: 22.09.2010 Сообщений: 18
	24.09.2010, 20:35 [ТС]
	Всем спасибо, я сам разобрался! 0