Форум программистов, компьютерный форум, киберфорум
Наши страницы
Алгоритмы
Войти
Регистрация
Восстановить пароль
 
Рейтинг 4.63/8: Рейтинг темы: голосов - 8, средняя оценка - 4.63
-Sky-
3 / 3 / 0
Регистрация: 22.09.2010
Сообщений: 18
1

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

22.09.2010, 19:08. Просмотров 1471. Ответов 10
Метки нет (Все метки)

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

Система:
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.
Помогите кто чем может
0
QA
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
22.09.2010, 19:08
Ответы с готовыми решениями:

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

какой способ лучше для базового изучения 8.2
Доброго времени суток! Я работаю с 7кой, по 8.1 есть некоторые теоретические знания. Посоветуйте,...

Какой способ измерения времени лучше использовать для получения данных с comport?
Нужно получать время между измерениями с comport что лучше, использовать? Точность нужна до 0,01с

Какой SSD подойдет для данной системы
Задача: подобрать оптимальный SDD для данной системы. Суть использования: OS, Photoshop,...

Ограничения для численного решения дифуров
Собственно вопрос, как задать ограничения для решения дифуров? Везде примеры, где условия вида...

10
KuKu
1560 / 1038 / 93
Регистрация: 17.04.2009
Сообщений: 2,995
22.09.2010, 21:23 2
аналитически попробуй , по размеру как раз выйдет как на один из методов, даже чуть меньше. ))
0
-Sky-
3 / 3 / 0
Регистрация: 22.09.2010
Сообщений: 18
22.09.2010, 21:33  [ТС] 3
Так в том то все и дело, что аналитически не получится! там слишком кривая ерунда получается...
0
KuKu
1560 / 1038 / 93
Регистрация: 17.04.2009
Сообщений: 2,995
22.09.2010, 22:33 4
Шаг попробуй поуменьшать. f1 и f2 у тебя какие то мутные еще...

Добавлено через 13 минут
это про рунге-кутта говорю, про остальные не знаю
0
22.09.2010, 22:33
-Sky-
3 / 3 / 0
Регистрация: 22.09.2010
Сообщений: 18
22.09.2010, 23:23  [ТС] 5
да ... там константы просто остались от старой системы, а так все нормально с ними
0
KuKu
1560 / 1038 / 93
Регистрация: 17.04.2009
Сообщений: 2,995
23.09.2010, 00:21 6
C++
1
result:=y;
вродь у рунге-кута решаем уравнение вида x'=f(x,t) и y=f(y,t) , что просто есть уравнения без дополнительных преобразований.
0
-Sky-
3 / 3 / 0
Регистрация: 22.09.2010
Сообщений: 18
23.09.2010, 00:39  [ТС] 7
Я уже говорил, что в функциях 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 секунд
то есть система которую мне нужно решить
0
KuKu
1560 / 1038 / 93
Регистрация: 17.04.2009
Сообщений: 2,995
23.09.2010, 03:38 8
метод рунге-кутта решает уравнение вида 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  [ТС] 9
ты ошибся, то что ты написал верно для одномерного случая, а у меня система из 2ух дифуров и у меня 3 переменных (пусть будут не t,x,y как у меня, а x,y,z), тогда в общем случае можно написать
dy/dx=f1(x,y,z)
dz/dx=f2(x,y,z)
тогда то и получается что ты написал
0
KuKu
1560 / 1038 / 93
Регистрация: 17.04.2009
Сообщений: 2,995
23.09.2010, 11:29 10
ну если так то извиняюсь. Ошибка хоть как будет при многом количестве итерраций. Опять же можно уменьшать шаг или попробуй описать приближения более высоких порядков. Можно попробывать часть решить аналитически, связь между x и y находится легко аналитчески это точно. Дальше решать только одно уравнение, ошибку это должно уменьшить, но не думаю что существенно. Можешь поставить математические пакеты и пусть они попробуют решить аналитически ... не надо пренебрегать теорие в данном случае, в аналитическом нет ошибок вообще. На форуме есть раздел дифуров, там думаю математиков найдешь, которые тебе скажут точно есть ли решение у подобных уравнений.
0
-Sky-
3 / 3 / 0
Регистрация: 22.09.2010
Сообщений: 18
24.09.2010, 20:35  [ТС] 11
Всем спасибо, я сам разобрался!
0
24.09.2010, 20:35
Answers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
24.09.2010, 20:35

Есть ли более простой и оптимальный способ решения данной задачи?
Добрый день! Ситуация следующая: Есть класс Event. public class Event { public string...

Какой цикл (For, While или Repeat) лучше использовать для решения задач с последовательностями?
Какой цикл For, While или Repeat лучше использовать для решения задач с последовательностями.

Какой способ лучше
Производятся вычисления в цикле, где в зависимости от условий подключаются или не подключаются...


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

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

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