State Vectorization for ODE 45

Hi all.

I'm trying to model a dyanmical system as following using vectorization of state for ODE45. My model includes three states that for two system there will be 6 states all in all. However, for some reason I'm gonna model these states by vectorization of states for solving by ODE45. However, as I checked the resutls, results are slightly different with respect to each other. Can you help why vectorization causes such these differences in the results? Thank you in Advance.

First code is as below:

clear all;clc;
%%
global b r I x0 
b=3.0;
r=0.02;
I=2.8;
x0=-1.6;
initial_condition=[-1,0.2,0.4,0.5,0,-0.2];
tspan=[0 3000];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t,y]=ode45(@EquationSys,tspan,initial_condition);
figure(1)
plot(t,y(:,1))
hold on
%%
function dy=EquationSys(t,y)
global b r I x0 
x1=y(1);
y1=y(2);
z1=y(3);
x2=y(4);
y2=y(5);
z2=y(6);
dy=[y1-x1^3+b*x1^2-z1+I;
    1-5*x1^2-y1;
    r*(4*(x1-x0)-z1);
    y2-x2^3+b*x2^2-z2+I;
    1-5*x2^2-y2;
    r*(4*(x2-x0)-z2);
    ];
end
%%

And Second code for vectorization is :

clear all;clc;
%%
global b r I x0 
b=3.0;
r=0.02;
I=2.8;
x0=-1.6;
initial_condition=[-1,0.2,0.4,0.5,0,-0.2];
tspan=[0 1000];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t,state]=ode45(@EquationSys,tspan,initial_condition);
figure(1)
plot(t,state(:,1))
hold on
%%
function dy=EquationSys(t,state)
global b r I x0 
x=state(1:2);
y=state(3:4);
z=state(5:6);
dy=[y-x.^3+b.*x.^2-z+I;
    1-5.*x.^2-y;
    r.*(4.*(x-x0)-z);
    ];
end
%%

[Shashi Kiran]

Hi @shahin sharafi,

After analyzing your code, I made some simple adjustments to the vectorized function to work as expected:

function dy=EquationSys(t,state)
global b r I x0 
x=state([1, 4]);
y=state([2, 5]);
z=state([3, 6]);
dy=zeros(6,1);
dy([1, 4]) = y - x.^3 + b.*x.^2 - z + I;
dy([2, 5]) = 1 - 5.*x.^2 - y;
dy([3, 6]) = r.*(4.*(x-x0) - z);
end

Hope this helps!

[shahin sharafi]

Thank you so much! it works

Answer 1

[Shivam Gothi]

Hello @shahin sharafi,

This happened because your state vector in case-1 is:

For the system of equations to give the same output their initial conditions should match. Therefore,just change the initial_condition vector in case 2 (Vactorised approach) as:

initial_condition=[-1,0.5,0.2,0,0.4,-0.2]; 

This will result in same solution for both the cases. I have attached the plot below. (Note:- there are two graphs, but they coincided exactly)

I hope this helps !

[Shivam Gothi]

Hello @shahin sharafi

sorry, I just made a typing mistake in the answer. Now it is corrected.