8 views (last 30 days)

Show older comments

I have a task to take 100 samples of a plot at a fixed value of parameters a2, omega1, and omega2 and I will be changing a1 and take 100 samples. It consumes so much time so I am wondering is their a way to take a 100 samples of a plot.

clear all; close all; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Simulation %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%avgnum=10;

%for avg=1:avgnum

OMEGA1= 1:0.2:6.4;

OMEGA2= 6.6:0.2:12;

%Erroravg_nonzero_x1_y1=zeros([1,length(OMEGA1)]);

%Erroravg_zero_x1_y1=zeros([1,length(OMEGA1)]);

%Erroravg_nonzero_x2_y2=zeros([1,length(OMEGA2)]);

%Erroravg_zero_x2_y2=zeros([1,length(OMEGA2)]);

for num1=1:length(OMEGA1)

for num2=1:length(OMEGA2)

%value of constants

a1=0.2;a2=0.3;

%omega1=5;omega2=4;

omega1=OMEGA1(num1);

omega2=OMEGA2(num2);

%omega1= 1:0.5:10;

%omega2= 11:0.5:20;

G=1;C12=0.01;C21=0.02;

dt=0.01; %step size

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1(1)=0.5;

y1(1)=0.5;

x2(1)=1;

y2(1)=1;

for i=2:1000

x1(i)=x1(i-1)+((a1-x1(i-1)^2-y1(i-1)^2)*x1(i-1)-omega1*y1(i-1)+G*C12*(x2(i-1)-x1(i-1)))*dt;

y1(i)=y1(i-1)+((a1-x1(i-1)^2-y1(i-1)^2)*y1(i-1)+omega1*x1(i-1)+G*C12*(y2(i-1)-y1(i-1)))*dt;

x2(i)=x2(i-1)+((a2-x2(i-1)^2-y2(i-1)^2)*x2(i-1)-omega2*y2(i-1)+G*C21*(x1(i-1)-x2(i-1)))*dt;

y2(i)=y2(i-1)+((a2-x2(i-1)^2-y2(i-1)^2)*y2(i-1)+omega2*x2(i-1)+G*C21*(y1(i-1)-y2(i-1)))*dt;

end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Observation %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ndatamin=4;Ndatamax=31;

Data=Ndatamin:1:Ndatamax;

Data=Data/Ndatamax;

Error_nonzero_x1=zeros([1,Ndatamax-Ndatamin+1]);%initialization

Error_zero_x1=zeros([1,Ndatamax-Ndatamin+1]);%initialization

Error_zero_y1=zeros([1,Ndatamax-Ndatamin+1]);%initialization

Error_nonzero_x2=zeros([1,Ndatamax-Ndatamin+1]);%initialization

Error_zero_x2=zeros([1,Ndatamax-Ndatamin+1]);%initialization

Error_nonzero_y2=zeros([1,Ndatamax-Ndatamin+1]);%initialization

Error_zero_y2=zeros([1,Ndatamax-Ndatamin+1]);%initialization

count=0;

for N_m=Ndatamin:Ndatamax

count=count+1;

N_measurements=N_m;

N_basis=31;

%index=randi([20000,30000],1,N_m);

index=randi([100,499],1,N_m);

Xdot1=zeros([1,N_measurements]);

Ydot1=zeros([1,N_measurements]);

Xdot2=zeros([1,N_measurements]);

Ydot2=zeros([1,N_measurements]);

for ni=1:N_measurements

Xdot1(ni)=(x1(index(ni)+1)-x1(index(ni)))/dt;

Ydot1(ni)=(y1(index(ni)+1)-y1(index(ni)))/dt;

Xdot2(ni)=(x2(index(ni)+1)-x2(index(ni)))/dt;

Ydot2(ni)=(y2(index(ni)+1)-y2(index(ni)))/dt;

end

M=zeros([N_measurements,N_basis]);

for i=1:N_measurements

for j=1:N_basis

if j==1

M(i,j)=1;

elseif j==2

M(i,j)=x1(index(i));

elseif j==3

M(i,j)=y1(index(i));

elseif j==4

M(i,j)=x2(index(i));

elseif j==5

M(i,j)=y2(index(i));

elseif j==6

M(i,j)=x1(index(i))^2;

elseif j==7

M(i,j)=x2(index(i))^2;

elseif j==8

M(i,j)=y1(index(i))^2;

elseif j==9

M(i,j)=y2(index(i))^2;

elseif j==10

M(i,j)=x1(index(i))*x2(index(i));

elseif j==11

M(i,j)=x1(index(i))*y1(index(i));

elseif j==12

M(i,j)=x1(index(i))*y2(index(i));

elseif j==13

M(i,j)=x2(index(i))*y1(index(i));

elseif j==14

M(i,j)=x2(index(i))*y2(index(i));

elseif j==15

M(i,j)=y1(index(i))*y2(index(i));

elseif j==16

M(i,j)=x1(index(i))^3;

elseif j==17

M(i,j)=y1(index(i))^3;

elseif j==18

M(i,j)=x2(index(i))^3;

elseif j==19

M(i,j)=y2(index(i))^3;

elseif j==20

M(i,j)=x1(index(i))^2*x2(index(i));

elseif j==21

M(i,j)=x1(index(i))^2*y1(index(i));

elseif j==22

M(i,j)=x1(index(i))^2*y2(index(i));

elseif j==23

M(i,j)=x2(index(i))^2*x1(index(i));

elseif j==24

M(i,j)=x2(index(i))^2*y1(index(i));

elseif j==25

M(i,j)=x2(index(i))^2*y2(index(i));

elseif j==26

M(i,j)=y1(index(i))^2*x1(index(i));

elseif j==27

M(i,j)=y1(index(i))^2*x2(index(i));

elseif j==28

M(i,j)=y1(index(i))^2*y2(index(i));

elseif j==29

M(i,j)=y2(index(i))^2*x1(index(i));

elseif j==30

M(i,j)=y2(index(i))^2*x2(index(i));

else

M(i,j)=y2(index(i))^2*y1(index(i));

end

end

end

end

figure

hold on

plot(x1,'r')

plot(x2,'g')

Sindar
on 30 Dec 2020

Edited: Sindar
on 30 Dec 2020

First, what are you actually doing with M, Xdot1, etc.? Currently, the code overwrites them each N_m loop, without doing anything with the results

Try this: define sampled vectors, so you can later vectorize things:

index=randi([100,499],1,N_m);

x1_samp0 = x1(index);

x1_samp1 = x1(index+1);

...

Xdot1=(x1_samp1-x1_samp0)/dt;

...

M(:,1)=1;

M(:,2)=x1_samp0;

...

M(:,6)=x1_samp0.^2;

edit: changed my mind on other option; don't do it. Here are some general tips, though:

- if you have a set of conditions that can be expressed as variable=first, variable=second, etc. -- use a switch instead of if-elseif
- If you have a loop running over a set of conditions and operations which have nothing else in common -- just run each line individually, then fill in the "else" portion after

for j=1:N_basis

if j==1

M(i,j)=1;

elseif j==2

M(i,j)=x1(index(i));

elseif j==3

M(i,j)=y1(index(i));

end

end

does the same as the below, but cleaner and faster:

M(i,1)=1;

M(i,2)=x1(index(i));

M(i,3)=y1(index(i));

- many, many, many for loops can be replaced with vector operations. It's kinda the point of Matlab

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

Start Hunting!