How do I solve differential equation in while/for loop using ODE45?

4 views (last 30 days)
Hi, I'm trying to simulate Rossler Oscillators with a 4 node system. . I would like to generalize this function by running this in a while/for loop, so that I can solve for 8-node, 12-node, n-node by just changing this in the function.
The initial conditions and ode function, etc. are in the main code. Below is the actual working code for simulating a 4 node system,
function yprime = rossler(t,y,K)
a = 0.18; b = 0.2; c = 5.7;
f = [
-(y(2)+y(3)); %dy1
(y(1) + a*y(2)); %dy2
b+y(3)*(y(1)-c); %dy3
-(y(5)+y(6)); %dy4
(y(4) + a*y(5)); %dy5
b+y(6)*(y(4)-c); %dy6
-(y(8)+y(9)); %dy7
(y(7) + a*y(8)); %dy8
b+y(9)*(y(7)-c); %dy9
-(y(11)+y(12)); %dy10
(y(10) + a*y(11)); %dy11
b+y(12)*(y(10)-c); %dy12
yprime = f - (K*y);
Below is the code I am trying to generalize by using for/while loop. Correct me if I am wrong. Thank you!
function yprime = cycfun(t,y,K)
a = 0.18; b = 0.2; c = 5.7;
n = 11; i = 1;
while i < n
i = i + 3;
f = [
-(y(i+1)+y(i+2)); %dy1
(y(i) + a*y(i+1)); %dy2
b+y(i+2)*(y(i)-c);]; %dy3
end
yprime = f - (K*y);
When I try running this code, I get matrix dimensions must agree. I know I'm messing it up somewhere, just trying to figure out! Thanks!
Below is the main code where I call the function:
clear all; clc;
tspan = [0,500];
options = odeset('stats', 'on');
%adjacency matrix, diagonal matrix, laplacian matrix
A = [0 1 1 1; 1 0 1 1; 1 1 0 0; 1 1 0 0];
D = [3 0 0 0; 0 3 0 0; 0 0 2 0; 0 0 0 2];
L = D - A;
%feedback structure, coupling factor, kronecker product
H = [1 0 0; 0 0 0; 0 0 0];
C = 0.001;
K = kron(L,H);
K = C*K;
%initial condition and ODE function
yzero = [0.093;0.3;0.065;0.093;0.3;0.065;0.093;0.3;0.065;0.093;0.3;0.065];
[T,Y] = ode45(@cycfun, tspan, yzero,options, K);
%plotting square difference
%--- NODE 1 & 2 ---%
%Square difference of node 1 and 2
q1 = ((Y(:,1) - Y(:,4)).^2) + ((Y(:,2) - Y(:,5)).^2) + ((Y(:,3) - Y(:,6)).^2);
w1 = q1;
%Plotting node 1 and 2
t1 = linspace(0,500,length(w1));
figure,plot(t1,w1);
title('Coupled Rossler Oscillator');
xlabel('t');
ylabel('(x11 - x21)^2 + (x12 - x22)^2 + (x13 - x23)^2');
%--- NODE 3 & 4 ---%
%Square difference of node 3 and 4
q2 = ((Y(:,7) - Y(:,10)).^2) + ((Y(:,8) - Y(:,11)).^2) + ((Y(:,9) - Y(:,12)).^2);
w2 = q2;
%Plotting node 1 and 2
t2 = linspace(0,500,length(w2));
figure,plot(t2,w2);
title('Coupled Rossler Oscillator');
xlabel('t');
ylabel('(x31 - x41)^2 + (x32 - x42)^2 + (x33 - x43)^2');
  4 Comments
KSSV
KSSV on 30 Nov 2017
I suspect problem is with your input...you should input y whose size should be 12...else error..show us how did you call the function?
Sreenath Umagandhi
Sreenath Umagandhi on 30 Nov 2017
@KSSV, thank you. I have updated my question with the code where I call the function.

Sign in to comment.

Accepted Answer

Torsten
Torsten on 30 Nov 2017
Edited: Torsten on 30 Nov 2017
function yprime = cycfun(t,y,K)
a = 0.18; b = 0.2; c = 5.7;
n = 12; i = 1; f = [];
while i < n
f = [f;-(y(i+1)+y(i+2));(y(i) + a*y(i+1)); b+y(i+2)*(y(i)-c)];
i = i+3;
end
yprime = f - (K*y);
Best wishes
Torsten.
  3 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!