清空筛选条件
清空筛选条件

Adjacency matrix for ODE45

42 次查看（过去 30 天）
N/A
N/A 2024-9-4，21:35
评论： N/A 2024-9-5，15:19
Hi
I am trying to model a systema as below:
where kij is :
I tried to model two sets of states manually first, without Aij (the adjacency matrix) as below:
clear all;clc;
%%
global b r I x0 alpha beta gamma
b=3.0;
r=0.006;
I=2.8;
x0=-1.6;
alpha=1;
beta=10;
gamma=0.5;
initial_condition=[-1,0.2,0.4,0.8,0,0.2,0.1,0.2];
tspan=[0 3000];
[t,y]=ode45(@EquationSys,tspan,initial_condition);
plot(t,y(:,1))
hold on
%%
function dy=EquationSys(t,y)
global b r I x0 alpha beta gamma
x1=y(1); %-1
y1=y(2);%0.2
z1=y(3);%0.4
x2=y(4);%0.8
y2=y(5);%0
z2=y(6);%0.2
k12=y(7);%0.1
k21=y(8);%0.2
dy=[y1-x1^3+b*x1^2-z1+I+k12*(x2-x1);
1-5*x1^2-y1;
r*(4*(x1-x0)-z1);
y2-x2^3+b*x2^2-z2+I+k21*(x1-x2);
1-5*x2^2-y2;
r*(4*(x2-x0)-z2);
k12*(alpha*exp(-beta*(x1 - x2)^2) - gamma*(k12+1));
k21*(alpha*exp(-beta*(x2 - x1)^2) - gamma*(k21+1));
];
end
%
%%%%%%%%%%%%%%%%%
For some reason, I want to add Aij to this model and vectorization the code and I programmed the code as below:
clear all;clc;
% HR Parameters
global b r x0 I N alpha beta gamma Adj_Matrix
b=3.0;
r=0.006;
I=2.8;
x0=-1.6;
alpha=1;
beta=10;
gamma=0.5;
%%
% Simulation parameters
T = 3000; % Total time
time_span = [0 T]; % Time vector
%%
% Number of neurons
% n=2; %Number of neurons on each edges
N = 2; % Number of all neurons
%% Adjacency Matrix
% Initialize the matrix with zeros
Adj_Matrix = zeros(N, N);
Adj_Matrix(1,2)=1;
Adj_Matrix(2,1)=1;
%%
% Flatten initial conditions and initial coupling strengths into a single vector
initial_conditions = [-1,0.8,0.2,0,0.4,0.2,0,0.1,0.2,0];
%%
% Solve the ODE for the coupled system with STDP
[time, sol] = ode45(@ hr_equations, time_span, initial_conditions);
%x, y, z, and k from the solution
x = sol(:, 1:N);
y = sol(:, N+1:2*N);
z = sol(:, 2*N+1:3*N);
k = sol(:, 3*N+1:end);
%%
% % Plotting results
%
% % If you want to visualize all the neurons' behavior together
plot(time,x(:,1))
hold on
%%
function dstate_dt = hr_equations(time, state)
global b r x0 I N alpha beta gamma Adj_Matrix
% Reshape the state vector
x = state(1:N);
y = state(N+1:2*N);
z = state(2*N+1:3*N);
k = reshape(state(3*N+1:end), N, N);
% HR Model Equations
dx = y-x.^3+b*x.^2-z+I+sum((Adj_Matrix.*k).*(x'-x),2);
dy = 1-5*x.^2-y;
dz = r*(4*(x-x0)-z);
dk = zeros(N, N);
for i = 1:N
for j = 1:N
if i ~= j
dk(i,j) = k(i,j) * (alpha * exp(-beta * (x(i) - x(j))^2) - gamma * (k(i,j) + 1));
end
end
end
% Flatten the derivatives to return as a single column vector
dstate_dt = [dx; dy; dz; reshape(dk, [], 1)];
end
However, after running this code and comparing two plots I found that there is slight differences between two results. I cannot find the issue. Can you help me in this regard?
Thanks

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-9-4，23:41
移动：Torsten 2024-9-4，23:42
Ran in:
k = [0,0.1,0.2,0];
reshape(k,2,2)
ans = 2x2
0 0.2000 0.1000 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
I think you want the reshaped matrix to be
reshape(k,2,2).'
ans = 2x2
0 0.1000 0.2000 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by