Решение системы дифуров. (Тема со звёздочкой)

@SpaceQuester · Регистрация: 19.01.2016

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

Добрый день.

Требуется решить систему ОДЕ. Классическое решение простое. Допустим у меня есть система из четырёх уравнений, пишем их в один файл с названием ODEsystem.m:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
function [ dy ] = ODEsystem(t, y)
 
k = 55;
 
dy = zeros(4, 1);
dy(1) = k * y(2) * y(3); % V
dy(2) = -y(1) * y(3); % m
dy(3) = -0.5 * y(1) * y(2); % n
dy(4) = -0.5 * y(1) * y(2); % h
 
end

В другом файле с названием ODEsystemSolver.m мы вызываем решалку и рисуем решение:

Matlab M
1
2
3
4
5
6
7
8
9
10
close all;
clc;
clear;
 
time = [0:0.1:2];
initial_states = [1 0 1 0];
 
[T, Y] = ode45(@ODEsystem, time, initial_states);
 
plot(T ,Y(:, 4), '-k');

Всё это работает.
Но есть одно большое НО. В данном случае четыре уравнения записаны в массив одной переменной dy, что очень плохо. Данный подход не работает когда требуется подсчитать взаимодействие нескольких систем, каждая из которых описывается системой нескольких уравнений. У меня допустим есть нода, (узел), она описывается четырьмя переменными: V, m, n, h. Когда нода одна, я могу это всё описать как dy(1), dy(2), dy(3), dy(4). Но когда нод несколько мне нужно записать каждую переменную отдельно с её индексом: V(i), m(i), n(i), h(i).

Вопрос такой: как это корректно прописать для MATLAB?

К сообщению приложен PDF документ, решать мы будем систему уравнений Ходжкина-Хаксли, где один нейрон описывается четырьмя уравнениями, но таких нейронов N штук, и связаны они разностной схемой: sum (V_j - V_i). Рассматривать будем на нижнем урезанном примере, что бы лучше разобраться в синтаксисе записи. ODEsolver.pdf

У кого есть мысли, как это прописать, пожалуйста, подскажите.

@SSC · 12.10.2017, 07:56

Запуская солвер ode45, Вы по сути сразу определяете число дифф. уравнений которые будете решать. Это число определяется по числу элементов в векторе начальных условий ( initial_states ).
Во время работы солвера изменить число уравнений нельзя, поэтому солвер запускается с максимальным количеством уравнений, или если необходимо изменить число решаемых уравнений останавливают этот запуск солвера, готовят данные для следующего этапа, и запускают солвер для продолжения решения.
По поводу функции ODEsystem, Вам никто не запрещает из нее обращаться к другим функциям и организовывать в ней циклические вычисления.
Основное условие, чтобы в заданный момент времени t, и известных величинах параметров у, Вы определили к моменту завершения функции величины всех производных dy.
Гибкость обработки может быть построена путем определения числа уравнений внутри функции следующим образом:

Matlab M
1
2
3
4
5
6
7
function [ dy ] = ODEsystem(t, y) 
k = 55;
n=length(y);
dy = zeros(n, 1);
for i=1:n
 dy(i)=.....
end

@SpaceQuester · 12.10.2017, 08:55 **[ТС]**

Спасибо за комментарий. Нет, на ходу я не буду менять число уравнений. Как делать динамически и задавать уравнения в цикле прямо в функции ODEsystem я разобрался, но такое решение покажу позже.

Пока текущие результаты показываю. Принудительно введены уравнения для двух нейронов.
Файл: ODEsystemSolver.m

Matlab M
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
close all;
clc;
clear;
 
time = [0:0.00001:1];
initial_states = [-61.5364 0.0789 0.3718 0.4723 34.3334 0.9165 0.5626 0.2440];
 
global I_app;
I_app_min = 5.20;
I_app_max = 5.34;
I_app = I_app_min + (I_app_max - I_app_min).*rand(2, 1);
 
tic;
initime = cputime;
time1   = clock;
 
[T, Y] = ode45(@ODEsystem, time, initial_states);
 
fintime = cputime;
elapsed = toc;
time2   = clock;
fprintf('TIC TOC: %g\n', elapsed);
fprintf('CPUTIME: %g\n', fintime - initime);
fprintf('CLOCK:   %g\n', etime(time2, time1));
 
figure('units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on');
plot(T ,Y(:, 1), 'LineWidth', 1);
hold on;
plot(T ,Y(:, 5), 'LineWidth', 1);
hold on;
grid on;
%xlim([0 10]);
%ylim([-80 40]);
set(gca, 'FontSize', 24);
title('V. MATLAB ode45 result', 'FontSize', 50);
xlabel('Time, sec', 'FontSize', 50);
ylabel('V, mV', 'FontSize', 50);
saveas(gcf, 'results_MATLAB.png');

Файл: ODEsystem.m

Matlab M
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
function [ dy ] = ODEsystem(t, y)
 
Node_count = 2;
Equations_per_node = 4;
Equations_count = Node_count * Equations_per_node;
 
C      = 1;
g_Na   = 120;
g_K    = 36;
g_Leak = 0.3;
E_Na   = 55;
E_K    = -77;
E_Leak = -54.4;
global I_app;
 
g_syn = 0.11;
k_syn = 0.2;
 
dy = zeros(Equations_count, 1);
 
dy(1) = 1000 * ( g_Na * y(2)^3 * y(4) * (E_Na - y(1)) + g_K * y(3)^4 * (E_K - y(1)) + g_Leak * (E_Leak - y(1)) + I_app(1) + g_syn * y(1) / (1 + exp(-y(5) / k_syn)) ) / C; % V
dy(2) = 1000 * ( (0.1 * (40 + y(1)) / (1 - exp(-(40+y(1))/10))) * (1 - y(2)) - (4 * exp(-(y(1)+65)/18)) * y(2) ); % m
dy(3) = 1000 * ( (0.01 * (55+y(1))/(1 - exp(-(55+y(1))/10))) * (1 - y(3)) - (0.125 * exp(-(y(1)+65)/80)) * y(3) ); % n
dy(4) = 1000 * ( (0.07 * exp(-(65+y(1))/20)) * (1 - y(4)) - (1/(exp(-(35+y(1))/10)+1)) * y(4) ); % h
 
dy(5) = 1000 * ( g_Na * y(6)^3 * y(8) * (E_Na - y(5)) + g_K * y(7)^4 * (E_K - y(5)) + g_Leak * (E_Leak - y(5)) + I_app(2) + g_syn * y(5) / (1 + exp(-y(1) / k_syn)) ) / C; % V
dy(6) = 1000 * ( (0.1 * (40 + y(5)) / (1 - exp(-(40+y(5))/10))) * (1 - y(6)) - (4 * exp(-(y(5)+65)/18)) * y(6)); % m
dy(7) = 1000 * ( (0.01 * (55+y(5))/(1 - exp(-(55+y(5))/10))) * (1 - y(7)) - (0.125 * exp(-(y(5)+65)/80)) * y(7)); % n
dy(8) = 1000 * ( (0.07 * exp(-(65+y(5))/20)) * (1 - y(8)) - (1/(exp(-(35+y(5))/10)+1)) * y(8)); % h
 
end

@SSC · 12.10.2017, 12:10

Вот вариант оформления Ваших уравнений

Matlab M
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
function [ dy ] = ODEsystem(t, y)
 
Node_count = 2;
Equations_per_node = 4;
Equations_count = Node_count * Equations_per_node;
 
C      = 1;
g_Na   = 120;
g_K    = 36;
g_Leak = 0.3;
E_Na   = 55;
E_K    = -77;
E_Leak = -54.4;
global I_app;
 
g_syn = 0.11;
k_syn = 0.2;
 
dy = zeros(Equations_count, 1);
 
YY=[y(5) y(1)];
 
for i=0:Equations_per_node:Equations_count-1
    dy(i+1) = 1000 * ( g_Na * y(i+2)^3 * y(i+4) * (E_Na - y(i+1)) + g_K * y(i+3)^4 * (E_K - y(i+1)) + g_Leak * (E_Leak - y(i+1)) + I_app(i/Equations_per_node+1) + g_syn * y(i+1) / (1 + exp(-YY(i/Equations_per_node+1) / k_syn)) ) / C; % V
    dy(i+2) = 1000 * ( (0.1 * (40 + y(i+1)) / (1 - exp(-(40+y(i+1))/10))) * (1 - y(i+2)) - (4 * exp(-(y(i+1)+65)/18)) * y(i+2) ); % m
    dy(i+3) = 1000 * ( (0.01 * (55+y(i+1))/(1 - exp(-(55+y(i+1))/10))) * (1 - y(i+3)) - (0.125 * exp(-(y(i+1)+65)/80)) * y(i+3) ); % n
    dy(i+4) = 1000 * ( (0.07 * exp(-(65+y(i+1))/20)) * (1 - y(i+4)) - (1/(exp(-(35+y(i+1))/10)+1)) * y(i+4) ); % h
end
end

@SpaceQuester · 12.10.2017, 12:48 **[ТС]**

Спасибо! Я только в отдельном for каждое уравнение генерировал. Пока вы не забыли тему про нейроны, нужно в первом уравнении хвост вот такой сделать sum_j^N ( g_syn * V(i) / (1 + exp(V(j)) / k_syn) ). Где V - это первые переменные в каждой подсистеме: y(0), y(5), y(9) и т.д.

Или лучше проще для понимания хвост сделать: g_syn * sum_j^N ( V(j) - V(i) )

@SSC · 12.10.2017, 13:09

У меня Ваших уравнений нет, поэтому то что вы пишете мне не понятно.
Прошлый раз я увидел одинаковость в группах по 4 уравнения, и отличия вынес как особенности в переменную YY

@SpaceQuester · 12.10.2017, 13:49 **[ТС]**

А теперь делаем еще круче. Для V_j делаем биологическое запаздывание в 7 миллисекунд. Т.е. для V_j берется значение из прошлого со сдвигом V_j(t - t_lag), где t_lag = 0.007.

Такая система дифуров уже не является однородной (ODE), а переходит в класс Delay Differential Equations (DDE). Что бы решить её воспользуемся функцией DDE в MATLAB:
https://www.mathworks.com/help... dde23.html
https://www.mathworks.com/help... elays.html

Теперь у нас три нейрона взаимодействуют друг с другом.
Файл с данными о задержками при t < 0: DDEhistory.m

Matlab M
1
2
function s = DDEhistory(t)
s = zeros(12, 1);

Файл с системой уравнений: DDEsystem.m

Matlab M
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
function [ dy ] = DDEsystem(t, y, Z)
 
Node_count = 3;
Equations_per_node = 4;
Equations_count = Node_count * Equations_per_node;
 
C      = 1;
g_Na   = 120;
g_K    = 36;
g_Leak = 0.3;
E_Na   = 55;
E_K    = -77;
E_Leak = -54.4;
global I_app;
 
g_syn = 0.15;
k_syn = 0.2;
 
dy = zeros(Equations_count, 1);
 
ylag1 = Z(:, 1);
ylag5 = Z(:, 2);
ylag9 = Z(:, 3);
 
dy(1) = 1000 * ( g_Na * y(2)^3 * y(4) * (E_Na - y(1)) + g_K * y(3)^4 * (E_K - y(1)) + g_Leak * (E_Leak - y(1)) + I_app(1) + g_syn * y(1) / (1 + exp(-ylag5(5) / k_syn)) + g_syn * y(1) / (1 + exp(-ylag9(9) / k_syn)) ) / C; % V
dy(2) = 1000 * ( (0.1 * (40 + y(1)) / (1 - exp(-(40+y(1))/10))) * (1 - y(2)) - (4 * exp(-(y(1)+65)/18)) * y(2) ); % m
dy(3) = 1000 * ( (0.01 * (55+y(1))/(1 - exp(-(55+y(1))/10))) * (1 - y(3)) - (0.125 * exp(-(y(1)+65)/80)) * y(3) ); % n
dy(4) = 1000 * ( (0.07 * exp(-(65+y(1))/20)) * (1 - y(4)) - (1/(exp(-(35+y(1))/10)+1)) * y(4) ); % h
 
dy(5) = 1000 * ( g_Na * y(6)^3 * y(8) * (E_Na - y(5)) + g_K * y(7)^4 * (E_K - y(5)) + g_Leak * (E_Leak - y(5)) + I_app(2) + g_syn * y(5) / (1 + exp(-ylag1(1) / k_syn)) + g_syn * y(5) / (1 + exp(-ylag9(9) / k_syn)) ) / C; % V
dy(6) = 1000 * ( (0.1 * (40 + y(5)) / (1 - exp(-(40+y(5))/10))) * (1 - y(6)) - (4 * exp(-(y(5)+65)/18)) * y(6)); % m
dy(7) = 1000 * ( (0.01 * (55+y(5))/(1 - exp(-(55+y(5))/10))) * (1 - y(7)) - (0.125 * exp(-(y(5)+65)/80)) * y(7)); % n
dy(8) = 1000 * ( (0.07 * exp(-(65+y(5))/20)) * (1 - y(8)) - (1/(exp(-(35+y(5))/10)+1)) * y(8)); % h
 
dy(9) = 1000 * ( g_Na * y(10)^3 * y(12) * (E_Na - y(9)) + g_K * y(11)^4 * (E_K - y(9)) + g_Leak * (E_Leak - y(9)) + I_app(3) + g_syn * y(9) / (1 + exp(-ylag1(1) / k_syn)) + g_syn * y(9) / (1 + exp(-ylag5(5) / k_syn)) ) / C; % V
dy(10) = 1000 * ( (0.1 * (40 + y(9)) / (1 - exp(-(40+y(9))/10))) * (1 - y(10)) - (4 * exp(-(y(9)+65)/18)) * y(10)); % m
dy(11) = 1000 * ( (0.01 * (55+y(9))/(1 - exp(-(55+y(9))/10))) * (1 - y(11)) - (0.125 * exp(-(y(9)+65)/80)) * y(11)); % n
dy(12) = 1000 * ( (0.07 * exp(-(65+y(9))/20)) * (1 - y(12)) - (1/(exp(-(35+y(9))/10)+1)) * y(12)); % h
 
end

Файл-решалка с отрисовкой результатов: DDEsystemSolver.m

Matlab M
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
close all;
clc;
clear;
 
time = [0:0.00001:5];
initial_states = [[-61.5364 0.0789 0.3718 0.4723] [-61.5364 0.0789 0.3718 0.4723] [34.3334 0.9165 0.5626 0.2440]];
 
global I_app;
I_app_min = 5.20;
I_app_max = 5.34;
I_app = I_app_min + (I_app_max - I_app_min).*rand(3, 1);
 
tic;
initime = cputime;
time1   = clock;
 
sol = dde23(@DDEsystem, [0.007, 0.007, 0.007] , @DDEhistory, time);
 
fintime = cputime;
elapsed = toc;
time2   = clock;
fprintf('TIC TOC: %g\n', elapsed);
fprintf('CPUTIME: %g\n', fintime - initime);
fprintf('CLOCK:   %g\n', etime(time2, time1));
 
figure('units', 'normalized', 'outerposition', [0 0 1 1], 'visible', 'on');
plot(sol.x ,sol.y(1, :), 'LineWidth', 1);
hold on;
plot(sol.x ,sol.y(5, :), 'LineWidth', 1);
hold on;
plot(sol.x ,sol.y(9, :), 'LineWidth', 1);
hold on;
grid on;
%xlim([0 10]);
%ylim([-80 40]);
set(gca, 'FontSize', 24);
title('V. MATLAB dde23 result', 'FontSize', 50);
xlabel('Time, sec', 'FontSize', 50);
ylabel('V, mV', 'FontSize', 50);
saveas(gcf, 'results_DDE_MATLAB.png');
 
clear all;

Что бы понять что за систему я программирую, приведу документ с описанием уравнений системы: ODEsolver.pdf

Новые блоги и статьи Все статьи Все блоги /
Учёным и волонтёрам проекта «Einstein@home» удалось обнаружить четыре гамма-лучевых пульсара в джете Млечного Пути Programma_Boinc 01.01.2026 Учёным и волонтёрам проекта «Einstein@home» удалось обнаружить четыре гамма-лучевых пульсара в джете Млечного Пути Сочетание глобально распределённой вычислительной мощности и инновационных. . .	Советы по крайней бережливости. Внимание, это ОЧЕНЬ длинный пост. 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? Ниже её машинный перевод. После долгих разбирательств я наконец-то вернула себе. . .	Thinkpad X220 Tablet — это лучший бюджетный ноутбук для учёбы, точка. Programma_Boinc 23.12.2025 Рецензия / Мнение/ Перевод Нашел на реддите интересную статью под названием The Thinkpad X220 Tablet is the best budget school laptop period . Ниже её машинный перевод. Thinkpad X220 Tablet —. . .	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/ нашёл схожую тему. Там приведён код на С++, который показывает только имена СОМ портов, типа,. . .

@SpaceQuester 2 / 2 / 0 Регистрация: 19.01.2016 Сообщений: 229
	12.10.2017, 12:48 [ТС]
	Спасибо! Я только в отдельном for каждое уравнение генерировал. Пока вы не забыли тему про нейроны, нужно в первом уравнении хвост вот такой сделать sum_j^N ( g_syn * V(i) / (1 + exp(V(j)) / k_syn) ). Где V - это первые переменные в каждой подсистеме: y(0), y(5), y(9) и т.д. Или лучше проще для понимания хвост сделать: g_syn * sum_j^N ( V(j) - V(i) ) 0

@SSC $Эксперт по математике/физике$ 3390 / 1913 / 571 Регистрация: 09.04.2015 Сообщений: 5,365
	12.10.2017, 13:09
	У меня Ваших уравнений нет, поэтому то что вы пишете мне не понятно. Прошлый раз я увидел одинаковость в группах по 4 уравнения, и отличия вынес как особенности в переменную YY 0