Converting MATLAB code into python code

21 次查看（过去 30 天）
aakash dewangan
aakash dewangan 2021-11-17
编辑： Mike Croucher 2021-11-21
Hello
I have written a MATLAB code, which is given below. I want to convert this code into pythod because I am not able to run this code on my laptop, which has limited Ram and very few Cores. So, i am going to ask high performance computation authority of my college to run this code. The problem is that they don't allow software which require GUI for running the codes. So, I want to convert this code into python to run it on clusters.
Please MATLAB community help me! MATLAB will always be my favourite :)
Thanks
tic
clc; close all; clear all;
syms p1(t) p2(t) p3(t) alpha_by_L beta_by_L gamma_by_L epsilon_by_L
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);
gamma_by_L_Matrix = [0.1:0.1:5];
beta_by_L_Matrix = [0.2:0.2:10];
alpha_by_L_Matrix = [0.1:0.1:5];
path = pwd ; % mention your path
myfolder = 'myfolder' ; % new folder name
folder = mkdir([path,filesep,myfolder]) ;
path = [path,filesep,myfolder] ;
parfor (iii = 1:1:length(gamma_by_L_Matrix))
gamma_by_L = gamma_by_L_Matrix(iii)
figure(1)
for jjj = 1:1:length(beta_by_L_Matrix)
beta_by_L = beta_by_L_Matrix(jjj)
for kkk = 1:1:length(alpha_by_L_Matrix)
alpha_by_L = alpha_by_L_Matrix(kkk)
% matrix terms
AA = 1/2 + (alpha_by_L)*(sin(pi*beta_by_L*t))^2;
BB = (alpha_by_L)*sin(pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);
CC = (alpha_by_L)*sin(pi*beta_by_L*t)*sin(3*pi*beta_by_L*t);
DD = 1/2 + (alpha_by_L)*(sin(2*pi*beta_by_L*t))^2;
EE = (alpha_by_L)*sin(3*pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);
FF = 1/2 + (alpha_by_L)*(sin(3*pi*beta_by_L*t))^2;
% matrix terms
GG = (pi^2)/2 + (gamma_by_L)*(sin(pi*beta_by_L*t))^2;
HH = (gamma_by_L)*sin(pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);
II = (gamma_by_L)*sin(3*pi*beta_by_L*t)*sin(pi*beta_by_L*t);
JJ = (4*pi^2)/2 + (gamma_by_L)*(sin(2*pi*beta_by_L*t))^2;
KK = (gamma_by_L)*sin(2*pi*beta_by_L*t)*sin(3*pi*beta_by_L*t);
LL = (9*pi^2)/2 + (gamma_by_L)*(sin(3*pi*beta_by_L*t))^2;
% RHS
MM = (epsilon_by_L)*sin(pi*beta_by_L*t);
NN = (epsilon_by_L)*sin(2*pi*beta_by_L*t);
OO = (epsilon_by_L)*sin(3*pi*beta_by_L*t);
Eq1 = AA*diff(p1,t,2) + BB*diff(p2,t,2) + CC*diff(p3,t,2) + GG*p1 + HH*p2 + II*p3 == 0; % Equation 1
Eq2 = BB*diff(p1,t,2) + DD*diff(p2,t,2) + EE*diff(p3,t,2) + HH*p1 + JJ*p2 + KK*p3 == 0; % Equation 2
Eq3 = CC*diff(p1,t,2) + EE*diff(p2,t,2) + FF*diff(p3,t,2) + II*p1 + KK*p2 + LL*p3 == 0; % Equation 3
[V,S] = odeToVectorField(Eq1, Eq2, Eq3);
ftotal = matlabFunction(V, 'Vars',{'t','Y'}); % Using readymade MATLAB function to solve using ODE 45
interval = [0 2/beta_by_L]; % Time Interval to run the program
for i = 1:6
if i == 1
y0 = [1 ;0 ;0 ;0 ;0 ;0 ] ; % First IC
elseif i == 2
y0 = [0 ;1 ;0 ;0 ;0 ;0 ] ; % Second IC
elseif i == 3
y0 = [0 ;0 ;1 ;0 ;0 ;0 ] ;
elseif i == 4
y0 = [0 ;0 ;0 ;1 ;0 ;0 ] ;
elseif i == 5
y0 = [0 ;0 ;0 ;0 ;1 ;0 ] ;
elseif i == 6
y0 = [0 ;0 ;0 ;0 ;0 ;1 ] ;
end
xSol = ode45( @(t,Y) ftotal(t,Y),interval,y0); % Using ODE 45 to solve state space form of ODE
tValues = linspace(interval(1),interval(2),1000); % dividing time interval
x1Values = deval(xSol,tValues,1); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
x2Values = deval(xSol,tValues,2); % number 2 denotes second solution likewise you can mention 2 ,3 & 4 for the next three solutions
x3Values = deval(xSol,tValues,3);
x4Values = deval(xSol,tValues,4);
x5Values = deval(xSol,tValues,5);
x6Values = deval(xSol,tValues,6);
if i == 1
col_11 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_12 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 2
col_21 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_22 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 3
col_31 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_32 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 4
col_41 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_42 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 5
col_51 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_52 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 6
col_61 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_62 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
end
end
A1 = [col_11 col_21 col_31 col_41 col_51 col_61];
A2 = [col_12 col_22 col_32 col_42 col_52 col_62] ;
A = (inv(A1))*A2 ;
eigval = eig(A);
abseigval = abs(eigval);
iteval = gamma_by_L
if (abseigval(1)<=1.00 && abseigval(2)<=1.00 && abseigval(3)<=1.00 && abseigval(4)<=1.00 && abseigval(5)<=1.00 && abseigval(6)<=1.00)
plot(beta_by_L,alpha_by_L,'k.', 'MarkerSize', 10)
xlabel('beta_by_L')
ylabel('alpha_by_L')
title("gamma_by_L " + gamma_by_L + " ")
hold on
else
plot(beta_by_L,alpha_by_L,'r.', 'MarkerSize', 10)
xlabel('beta_by_L')
ylabel('alpha_by_L')
title("gamma_by_L " + gamma_by_L + " ")
hold on
end
end
end
hold off
xlim([0 10])
ylim([0 5])
iii
end
toc
  2 个评论
Mike Croucher
Mike Croucher 2021-11-19
You don't need a GUI to run MATLAB code. You can run it in batch mode just like you can Python. A huge part of my job is to work with researchers who want to run MATLAB effectively on High Performance Computing clusters.
Your first issue with this code is 'What results do you want to save?' I'm guessing the variables beta_by_L,alpha_by_L since you are plotting them?
aakash dewangan
aakash dewangan 2021-11-19
Thank you very much Mike for reply :)
I am glad to know that I can use MATLAB in high performance computation clusters :)
From this code, I will be getting plots/graphs of "beta_by_L vs alpha_by_L" for different values of "gamma_by_L". For different values of gamma_by_L, I will get different graphs. I want all those graphs in my result.
Please help me and guide me to run this on high performance computing cluster.
I really need your help.

请先登录，再进行评论。

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

采纳的回答

Mike Croucher
Mike Croucher 2021-11-20
编辑：Mike Croucher 2021-11-21
Hi Aakash
To get this running on a HPC system, it would be beneficial to know the details of the HPC system. As I said in my comment to your question, you can run MATLAB in batch mode (that is without a GUI). The relevant section of the documentation is Start MATLAB program from Linux system prompt - MATLAB - MathWorks United Kingdom and the command switch you want is matlab -batch yourcode.m
Before going further down that route, however, I suggest some changes to your code.
Note that I have only done limited testing here and did the smallest amount of changes I could to get the speed-ups. You should check that it gives the results you need!
There's more that could be done to improve the readability of this code that might not affect the speed. I deliberately didn't do any of that.
On my machine, the first iteration of iii went down from 2942 to 20 seconds. A speed-up of 147 times.
The full run took 227 seconds to produce 50 plots on an 8 core system.
The only 'test' I've done is to check that the first plot given corresponding to iii=1 is the same in both your original code and the new one....and this was done by eye. Please check yourself that you get the results you expect.
The original code gave the following in 2942 seconds for iii=1
the final version of the code gave the following in 20 seconds for iii=1
The first idea is to move all of the slow symbolic calculations outside of the main loop. This provides most of the speed-up, around 70x or so
%% new version of code that speeds up the symbolic section
%% this is not the fastest version so far
tic
clc; close all; clear all;
syms p1(t) p2(t) p3(t) alpha_by_L beta_by_L gamma_by_L epsilon_by_L
syms Sgamma_by_L Sbeta_by_L Salpha_by_L
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);
gamma_by_L_Matrix = [0.1:0.1:5];
beta_by_L_Matrix = [0.2:0.2:10];
alpha_by_L_Matrix = [0.1:0.1:5];
%%You created an output folder here but you are not saving
%%anything to it. I've chaged it a little and save the figure of each
%%iteration.
resultsfolder = "myfolder";
mkdir(resultsfolder);
% Create the symbolic matrix, and a slightly modified ftotal once instead
% of every iteration. The new ftotal function can be used every iteration
% without being created from scratch
diff_p1_t_2 = diff(p1,t,2);
diff_p2_t_2 = diff(p2,t,2);
diff_p3_t_2 = diff(p3,t,2);
AA = 1/2 + (Salpha_by_L)*(sin(pi*Sbeta_by_L*t))^2;
BB = (Salpha_by_L)*sin(pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
CC = (Salpha_by_L)*sin(pi*Sbeta_by_L*t)*sin(3*pi*Sbeta_by_L*t);
DD = 1/2 + (Salpha_by_L)*(sin(2*pi*Sbeta_by_L*t))^2;
EE = (Salpha_by_L)*sin(3*pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
FF = 1/2 + (Salpha_by_L)*(sin(3*pi*Sbeta_by_L*t))^2;
% matrix terms
GG = (pi^2)/2 + (Sgamma_by_L)*(sin(pi*Sbeta_by_L*t))^2;
HH = (Sgamma_by_L)*sin(pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
II = (Sgamma_by_L)*sin(3*pi*Sbeta_by_L*t)*sin(pi*Sbeta_by_L*t);
JJ = (4*pi^2)/2 + (Sgamma_by_L)*(sin(2*pi*Sbeta_by_L*t))^2;
KK = (Sgamma_by_L)*sin(2*pi*Sbeta_by_L*t)*sin(3*pi*Sbeta_by_L*t);
LL = (9*pi^2)/2 + (Sgamma_by_L)*(sin(3*pi*Sbeta_by_L*t))^2;
% RHS
MM = (epsilon_by_L)*sin(pi*Sbeta_by_L*t);
NN = (epsilon_by_L)*sin(2*pi*Sbeta_by_L*t);
OO = (epsilon_by_L)*sin(3*pi*Sbeta_by_L*t);
Eq1 = AA*diff_p1_t_2 + BB*diff_p2_t_2 + CC*diff_p3_t_2 + GG*p1 + HH*p2 + II*p3 == 0; % Equation 1
Eq2 = BB*diff_p1_t_2 + DD*diff_p2_t_2 + EE*diff_p3_t_2 + HH*p1 + JJ*p2 + KK*p3 == 0; % Equation 2
Eq3 = CC*diff_p1_t_2 + EE*diff_p2_t_2 + FF*diff_p3_t_2 + II*p1 + KK*p2 + LL*p3 == 0; % Equation 3
[V,S] = odeToVectorField(Eq1, Eq2, Eq3);
ftotal = matlabFunction(V, 'Vars',{'t','Y','Salpha_by_L','Sbeta_by_L','Sgamma_by_L'}); % Using readymade MATLAB function to solve using ODE 45
parfor iii = 1:length(gamma_by_L_Matrix)
gamma_by_L = gamma_by_L_Matrix(iii);
myfig = figure(1)
for jjj = 1:1:length(beta_by_L_Matrix)
beta_by_L = beta_by_L_Matrix(jjj);
for kkk = 1:1:length(alpha_by_L_Matrix)
alpha_by_L = alpha_by_L_Matrix(kkk);
interval = [0 2/beta_by_L]; % Time Interval to run the program
%% This section could be rewritten a lot but I leave it the same
% to make it easier for you to follow the changes done so far
for i = 1:6
if i == 1
y0 = [1 ;0 ;0 ;0 ;0 ;0 ]; % First IC
elseif i == 2
y0 = [0 ;1 ;0 ;0 ;0 ;0 ]; % Second IC
elseif i == 3
y0 = [0 ;0 ;1 ;0 ;0 ;0 ];
elseif i == 4
y0 = [0 ;0 ;0 ;1 ;0 ;0 ];
elseif i == 5
y0 = [0 ;0 ;0 ;0 ;1 ;0 ];
elseif i == 6
y0 = [0 ;0 ;0 ;0 ;0 ;1 ];
end
xSol = ode45( @(t,Y) ftotal(t,Y,alpha_by_L,beta_by_L,gamma_by_L),interval,y0); % Using ODE 45 to solve state space form of ODE
tValues = linspace(interval(1),interval(2),1000); % dividing time interval
x1Values = deval(xSol,tValues,1); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
x2Values = deval(xSol,tValues,2); % number 2 denotes second solution likewise you can mention 2 ,3 & 4 for the next three solutions
x3Values = deval(xSol,tValues,3);
x4Values = deval(xSol,tValues,4);
x5Values = deval(xSol,tValues,5);
x6Values = deval(xSol,tValues,6);
if i == 1
col_11 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_12 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 2
col_21 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_22 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 3
col_31 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_32 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 4
col_41 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_42 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 5
col_51 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_52 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 6
col_61 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_62 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
end
end
A1 = [col_11 col_21 col_31 col_41 col_51 col_61];
A2 = [col_12 col_22 col_32 col_42 col_52 col_62] ;
A = (inv(A1))*A2 ;
eigval = eig(A);
abseigval = abs(eigval);
iteval = gamma_by_L;
%% You make a plot here for every iteration of the iii loop but are not saving it.
%% Is it these plots that you want to save?
if (abseigval(1)<=1.00 && abseigval(2)<=1.00 && abseigval(3)<=1.00 && abseigval(4)<=1.00 && abseigval(5)<=1.00 && abseigval(6)<=1.00)
plot(beta_by_L,alpha_by_L,'k.', 'MarkerSize', 10)
xlabel('beta_by_L')
ylabel('alpha_by_L')
title("gamma_by_L " + gamma_by_L + " ")
hold on
else
plot(beta_by_L,alpha_by_L,'r.', 'MarkerSize', 10)
xlabel('beta_by_L')
ylabel('alpha_by_L')
title("gamma_by_L " + gamma_by_L + " ")
hold on
end
end
end
hold off
%% An example of how to save the plots to a unique filename per iteration
filename = resultsfolder+filesep+"result_"+iii+".fig";
savefig(filename);
%% I would suggest saving a numerical result as well per iteration. Have a think about what that should be
xlim([0 10])
ylim([0 5])
iii
end
toc
The next thing you can do is vectorise your calls to deval. i.e. change
% original deval section
x1Values = deval(xSol,tValues,1); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
x2Values = deval(xSol,tValues,2); % number 2 denotes second solution likewise you can mention 2 ,3 & 4 for the next three solutions
x3Values = deval(xSol,tValues,3);
x4Values = deval(xSol,tValues,4);
x5Values = deval(xSol,tValues,5);
x6Values = deval(xSol,tValues,6);
if i == 1
col_11 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_12 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 2
col_21 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_22 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 3
col_31 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_32 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 4
col_41 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_42 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 5
col_51 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_52 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 6
col_61 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_62 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
end
to
%% New deval section
xvalues = deval(xSol,tValues,1:6);
if i == 1
col_11 = xvalues(:,1);
col_12 = xvalues(:,end);
elseif i == 2
col_21 = xvalues(:,1);
col_22 = xvalues(:,end);
elseif i == 3
col_31 = xvalues(:,1);
col_32 = xvalues(:,end);
elseif i == 4
col_41 = xvalues(:,1);
col_42 = xvalues(:,end);
elseif i == 5
col_51 = xvalues(:,1);
col_52 =xvalues(:,end);
elseif i == 6
col_61 = xvalues(:,1);
col_62 = xvalues(:,end);
end
This gave another 2x speed-up on my machine. One iteration of iii took 20 seconds now with the full calculation of 50 iterations taking 227 seconds on my 8 core machine. The full version of the new code is:
%% The fastest version so far
tic
clc; close all; clear all;
syms p1(t) p2(t) p3(t) alpha_by_L beta_by_L gamma_by_L epsilon_by_L
syms Sgamma_by_L Sbeta_by_L Salpha_by_L
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);
gamma_by_L_Matrix = [0.1:0.1:5];
beta_by_L_Matrix = [0.2:0.2:10];
alpha_by_L_Matrix = [0.1:0.1:5];
%%You created an output folder here but you are not saving
%%anything to it. I've chaged it a little and save the figure of each
%%iteration.
resultsfolder = "myfolder";
mkdir(resultsfolder);
% Create the symbolic matrix, and a slightly modified ftotal once instead
% of every iteration. The new ftotal function can be used every iteration
% without being created from scratch
diff_p1_t_2 = diff(p1,t,2);
diff_p2_t_2 = diff(p2,t,2);
diff_p3_t_2 = diff(p3,t,2);
AA = 1/2 + (Salpha_by_L)*(sin(pi*Sbeta_by_L*t))^2;
BB = (Salpha_by_L)*sin(pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
CC = (Salpha_by_L)*sin(pi*Sbeta_by_L*t)*sin(3*pi*Sbeta_by_L*t);
DD = 1/2 + (Salpha_by_L)*(sin(2*pi*Sbeta_by_L*t))^2;
EE = (Salpha_by_L)*sin(3*pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
FF = 1/2 + (Salpha_by_L)*(sin(3*pi*Sbeta_by_L*t))^2;
% matrix terms
GG = (pi^2)/2 + (Sgamma_by_L)*(sin(pi*Sbeta_by_L*t))^2;
HH = (Sgamma_by_L)*sin(pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
II = (Sgamma_by_L)*sin(3*pi*Sbeta_by_L*t)*sin(pi*Sbeta_by_L*t);
JJ = (4*pi^2)/2 + (Sgamma_by_L)*(sin(2*pi*Sbeta_by_L*t))^2;
KK = (Sgamma_by_L)*sin(2*pi*Sbeta_by_L*t)*sin(3*pi*Sbeta_by_L*t);
LL = (9*pi^2)/2 + (Sgamma_by_L)*(sin(3*pi*Sbeta_by_L*t))^2;
% RHS
MM = (epsilon_by_L)*sin(pi*Sbeta_by_L*t);
NN = (epsilon_by_L)*sin(2*pi*Sbeta_by_L*t);
OO = (epsilon_by_L)*sin(3*pi*Sbeta_by_L*t);
Eq1 = AA*diff_p1_t_2 + BB*diff_p2_t_2 + CC*diff_p3_t_2 + GG*p1 + HH*p2 + II*p3 == 0; % Equation 1
Eq2 = BB*diff_p1_t_2 + DD*diff_p2_t_2 + EE*diff_p3_t_2 + HH*p1 + JJ*p2 + KK*p3 == 0; % Equation 2
Eq3 = CC*diff_p1_t_2 + EE*diff_p2_t_2 + FF*diff_p3_t_2 + II*p1 + KK*p2 + LL*p3 == 0; % Equation 3
[V,S] = odeToVectorField(Eq1, Eq2, Eq3);
ftotal = matlabFunction(V, 'Vars',{'t','Y','Salpha_by_L','Sbeta_by_L','Sgamma_by_L'}); % Using readymade MATLAB function to solve using ODE 45
parfor iii = 1:length(gamma_by_L_Matrix)
gamma_by_L = gamma_by_L_Matrix(iii);
myfig = figure(1)
for jjj = 1:1:length(beta_by_L_Matrix)
beta_by_L = beta_by_L_Matrix(jjj);
for kkk = 1:1:length(alpha_by_L_Matrix)
alpha_by_L = alpha_by_L_Matrix(kkk);
interval = [0 2/beta_by_L]; % Time Interval to run the program
%% This section could be rewritten a lot but I leave it the same to make it easier to follow the changes done so far
for i = 1:6
if i == 1
y0 = [1 ;0 ;0 ;0 ;0 ;0 ]; % First IC
elseif i == 2
y0 = [0 ;1 ;0 ;0 ;0 ;0 ]; % Second IC
elseif i == 3
y0 = [0 ;0 ;1 ;0 ;0 ;0 ];
elseif i == 4
y0 = [0 ;0 ;0 ;1 ;0 ;0 ];
elseif i == 5
y0 = [0 ;0 ;0 ;0 ;1 ;0 ];
elseif i == 6
y0 = [0 ;0 ;0 ;0 ;0 ;1 ];
end
xSol = ode45( @(t,Y) ftotal(t,Y,alpha_by_L,beta_by_L,gamma_by_L),interval,y0); % Using ODE 45 to solve state space form of ODE
tValues = linspace(interval(1),interval(2),1000); % dividing time interval
xvalues = deval(xSol,tValues,1:6);
if i == 1
col_11 = xvalues(:,1);
col_12 = xvalues(:,end);
elseif i == 2
col_21 = xvalues(:,1);
col_22 = xvalues(:,end);
elseif i == 3
col_31 = xvalues(:,1);
col_32 = xvalues(:,end);
elseif i == 4
col_41 = xvalues(:,1);
col_42 = xvalues(:,end);
elseif i == 5
col_51 = xvalues(:,1);
col_52 =xvalues(:,end);
elseif i == 6
col_61 = xvalues(:,1);
col_62 = xvalues(:,end);
end
end
A1 = [col_11 col_21 col_31 col_41 col_51 col_61];
A2 = [col_12 col_22 col_32 col_42 col_52 col_62] ;
A = (inv(A1))*A2 ;
eigval = eig(A);
abseigval = abs(eigval);
iteval = gamma_by_L;
%% You make a plot here for every iteration of the iii loop but are not saving it.
%% Is it these plots that you want to save?
if (abseigval(1)<=1.00 && abseigval(2)<=1.00 && abseigval(3)<=1.00 && abseigval(4)<=1.00 && abseigval(5)<=1.00 && abseigval(6)<=1.00)
plot(beta_by_L,alpha_by_L,'k.', 'MarkerSize', 10)
xlabel('beta_by_L')
ylabel('alpha_by_L')
title("gamma_by_L " + gamma_by_L + " ")
hold on
else
plot(beta_by_L,alpha_by_L,'r.', 'MarkerSize', 10)
xlabel('beta_by_L')
ylabel('alpha_by_L')
title("gamma_by_L " + gamma_by_L + " ")
hold on
end
end
end
hold off
%% An example of how to save the plots to a unique filename per iteration
filename = resultsfolder+filesep+"result_"+iii+".fig";
savefig(filename);
%% I would suggest saving a numerical result as well per iteration. Have a think about what that should be
xlim([0 10])
ylim([0 5])
iii
end
toc
Hope this helps,
Mike
  1 个评论
aakash dewangan
aakash dewangan 2021-11-20
O My GOD!!!! O My GOD!!!! O My GOD!!!! O My GOD!!!! O My GOD!!!! O My GOD!!!! O My GOD!!!!
I can't express my feeling in words. You are legend in the MATLAB sir. Thank you so much for your help. thank you very very very very very much sir. It works exactly in the same way, and the results are also same. I don't know how to thank you for this help. You have no idea what you have done for me. MATLAB is best...... And you are the best in solving the problems of MATLAB.
Best Regards,

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by