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

Question

Sreenath Umagandhi 2017-11-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/370113-how-do-i-solve-differential-equation-in-while-for-loop-using-ode45

评论： Sreenath Umagandhi 2017-11-30

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 个评论
显示 2更早的评论隐藏 2更早的评论

KSSV 2017-11-30

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 2017-11-30

@KSSV, thank you. I have updated my question with the code where I call the function.

请先登录，再进行评论。

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

Answer 1

Torsten 2017-11-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/370113-how-do-i-solve-differential-equation-in-while-for-loop-using-ode45#answer_293918

编辑：Torsten 2017-11-30

在 MATLAB Online 中打开

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);