С Новым годом! Форум программистов, компьютерный форум, киберфорум
Delphi для начинающих
Форум программистов и сисадминов Киберфорум > Форум программистов > Delphi > Delphi для начинающих
Войти
Регистрация
Восстановить пароль
Блоги Сообщество Поиск Заказать работу  
 
Рейтинг 4.60/5: Рейтинг темы: голосов - 5, средняя оценка - 4.60
 Аватар для -Sky-
3 / 3 / 0
Регистрация: 22.09.2010
Сообщений: 18

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

22.09.2010, 18:55. Показов 1054. Ответов 1
Метки нет (Все метки)
Студворк — интернет-сервис помощи студентам
Привет!)
У меня проблемка, не могу решить системы дифуров численно.

Система:
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
IT_Exp
Эксперт
34794 / 4073 / 2104
Регистрация: 17.06.2006
Сообщений: 32,602
Блог
 22.09.2010, 18:55
Ответы с готовыми решениями:

Какой способ лучше для решения данной системы дифуров?
Привет!) У меня проблемка, не могу решить системы дифуров численно. Система: dx/dt=a1*x-b1*x*y dy/dt=-a2*x+b2*x*y где...

Ограничения для численного решения дифуров
Собственно вопрос, как задать ограничения для решения дифуров? Везде примеры, где условия вида y(0)=0 и \dot{y}(0) = 0, а в моем случае...

Какой метод решения у заданий?
1)Создать графическое приложение Frame содержащее меню выбора текстового файла, поля для ввода текста и кнопку Apply. Реализовать...

1
215 / 215 / 20
Регистрация: 18.05.2010
Сообщений: 865
22.09.2010, 20:41
Вы задеваете вопросы которые в большей части относятся к математике, эти вопросы задайте в разделе математики они вам помогут.
1
Надоела реклама? Зарегистрируйтесь и она исчезнет полностью.
BasicMan
Эксперт
29316 / 5623 / 2384
Регистрация: 17.02.2009
Сообщений: 30,364
Блог
 22.09.2010, 20:41
Помогаю со студенческими работами здесь

Решение системы дифуров.
По электротехнике задали решить систему уравнений в MATLAB'е, а я в нем практически не шарю. Помогите кто знает как это сделать. Дана...

Матричный метод решения системы
Матричным методом решить системы.Что не так?? Не могу понять как правельно сделать. k=4; U=; L=; ??? Error using ==> inv...

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

Определить метод решения системы уравнений
Здравствуйте! Не могу определить метод решения системы уравнений. По идее, это же должно решаться методом given... find?

Рекурсивный метод решения системы уравнений
Решить систему уравнений рекурсивным методом: a,b,q - вводятся пользователем


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

Или воспользуйтесь поиском по форуму:
2
Ответ Создать тему
Новые блоги и статьи
Модель микоризы: классовый агентный подход 3
anaschu 06.01.2026
aa0a7f55b50dd51c5ec569d2d10c54f6/ O1rJuneU_ls https:/ / vkvideo. ru/ video-115721503_456239114
 Owen Logic: О недопустимости использования связки «аналоговый ПИД» + RegKZR
ФедосеевПавел 06.01.2026
Owen Logic: О недопустимости использования связки «аналоговый ПИД» + RegKZR ВВЕДЕНИЕ Введу сокращения: аналоговый ПИД — ПИД регулятор с управляющим выходом в виде числа в диапазоне от 0% до. . .
Модель микоризы: классовый агентный подход 2
anaschu 06.01.2026
репозиторий https:/ / github. com/ shumilovas/ fungi ветка по-частям. коммит Create переделка под биомассу. txt вход sc, но sm считается внутри мицелия. кстати, обьем тоже должен там считаться. . . .
Расчёт токов в цепи постоянного тока
igorrr37 05.01.2026
/ * Дана цепь постоянного тока с сопротивлениями и напряжениями. Надо найти токи в ветвях. Программа составляет систему уравнений по 1 и 2 законам Кирхгофа и решает её. Последовательность действий:. . .
Новый CodeBlocs. Версия 25.03
palva 04.01.2026
Оказывается, недавно вышла новая версия CodeBlocks за номером 25. 03. Когда-то давно я возился с только что вышедшей тогда версией 20. 03. С тех пор я давно снёс всё с компьютера и забыл. Теперь. . .
Модель микоризы: классовый агентный подход
anaschu 02.01.2026
Раньше это было два гриба и бактерия. Теперь три гриба, растение. И на уровне агентов добавится между грибами или бактериями взаимодействий. До того я пробовал подход через многомерные массивы,. . .
Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост.
Programma_Boinc 28.12.2025
Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. Налог на собак: https:/ / **********/ gallery/ V06K53e Финансовый отчет в Excel: https:/ / **********/ gallery/ bKBkQFf Пост отсюда. . .
Кто-нибудь знает, где можно бесплатно получить настольный компьютер или ноутбук? США.
Programma_Boinc 26.12.2025
Нашел на реддите интересную статью под названием Anyone know where to get a free Desktop or Laptop? Ниже её машинный перевод. После долгих разбирательств я наконец-то вернула себе. . .
КиберФорум - форум программистов, компьютерный форум, программирование
Поддержать - Реклама - Условия использования - Обратная связь
Powered by vBulletin
Copyright ©2000 - 2026, CyberForum.ru