Calling a function iteratively within a for loop with various size outputs

2 visualizzazioni (ultimi 30 giorni)
Hi this could be a long one, thanks for your consdieration. I have a "for" loop which iterates over a growing differential equation (the inputs of thuis may or may not be known a-priori). I then iterate over each extension of the differential equation and solve it. But I write the differential equation to a .m file which I then call as a function and solve. The output of this differential equation is a matrix N rows and M columns, The number of rows (1501) shouldn't change between each iteration however the number of columns M will. In this example it should grow in steps of 3, 6, 9, 12, 15. This could change later though so I wouldn't like them constrained. For some reason the code goes through the first iteration then at the line:
A=CME_solver(maxBeta);
The code breaks.
To explain the script I first decide what order of extension with the for loop. I then provide some logic where I select certain combinations of a vector.
I then add and subtract these elements of a vector to generate a matrix of sums and differences of each combination based on where each element came from in the original vector.
By comparing the numbers in the matrix with the numbers in the initial vector I can extract all the numbers (and locations) from the matrix for all numbes which existed in the initial vector and thus match what sum and difference combinations produced them.
I have to then convert this to a differential equation which is written as a string in two .m files with the names CME_solver.m and CME_solver_no_coeff.m. These Are standard ode45's essentially written programaticaly and then placed into the same folder as the script.
I then want to get the result of this differential equation and store each result for each iteration in the cell array I and I_no_coeff which i would finally pass back from this function to a main script once the for loop is finished.
Please see code below (uncomment the first inputs here to run as a stand alone script if necesary)
function [I,I_no_coeff,f_CME] = CMEj_solver_writer(param,N,maxBeta,Ap0,As0,up_folder)
% up_folder=/home/[your_folder_here]/..;
% Ap0=0.9502;
% As0=0.0801;
% param(5)=12;
% param(7)=7.2;
% N=1500;
% maxBeta=0.433;
CME_extension = 5;
for harmonic=1:CME_extension
wp = 2*pi*param(5)*10^9;
ws = 2*pi*param(7)*10^9;
wi = wp-ws;
% wdel = abs(ws-wi)
w = [wp,ws,wi];
w=sort(w)
wc=harmonic*wp;
for k=1:harmonic^2
for j = 1:length(w)
for i=1:length(w);
a = w(j) + w(i);
b = w(j) - w(i);
check1 = ismembertol(a,w,0.001);
check2 = ismembertol(b,w,0.001);
if mod(w(j),wp)==0 || mod(w(i),wp)==0 && b~=0
if check1==0 && a<wc+1e6
w = [w,a];
end
if check2==0 && b>0
w = [b,w];
end
end
end
end
end
%toDelete=w>wc;
%w(toDelete)=[];
w=sort(w);
f=w./(2*pi)
fa=a./(2*pi)
for j=1:length(w)
for i=1:length(w)
sfg = w(j) + w(i);
check = ismembertol(sfg,w,0.01);
if check==1;
sfg_matrix(j,i) = sfg;
else
sfg_matrix(j,i) = 0;
end
if i>j
sfg_matrix(j,i) = 0;
end
dfg = w(j) - w(i);
check = ismembertol(dfg,w,0.01);
if check==1;
dfg_matrix(j,i) = dfg;
else
dfg_matrix(j,i) = 0;
end
end
end
dfg_trans = dfg_matrix';
mixing_matrix = sfg_matrix + dfg_trans;
for i = 1:length(w)
w_string{i} = ['w',num2str(i)];
k_string{i} = ['k',num2str(i)];
amp_string{i} = ['A(',num2str(i),')'];
end
for j=1:length(w)
for i=1:length(w)
%%%% mixing_matrix(j,i) for j<i is dfg ie top right half
%%%% mixing_matrix(j,i) for j<i is sfg ie bottom left half
%%%% mixing_matrix(j,i) for j=i is SHG of ji = jj = ii is NW-SE
%%%% diagonal
if ismembertol(mixing_matrix(j,i),w,0.001)==1
if j<i
[check,place] = ismembertol(mixing_matrix(j,i),w);
dA_string{place,i} = ['+ (',k_string{i},'-',k_string{j},')*',k_string{i},'*',k_string{j},'*',amp_string{i},'*conj(',amp_string{j},')*exp(1i*(',k_string{i},'-',k_string{j},')*x) '];
end
if j>i
[check,place] = ismembertol(mixing_matrix(j,i),w,0.001);
dA_string{place,i} = ['- (',k_string{i},'+',k_string{j},')*',k_string{i},'*',k_string{j},'*',amp_string{i},'*',amp_string{j},'.*exp(1i*(',k_string{i},'+',k_string{j},')*x) '];
end
if j==i
if ismembertol(mixing_matrix(j,i),w,0.001)==1
[check,place] = ismembertol(mixing_matrix(j,i),w) ;
dA_string{i+j,i+j} = ['- 0.5*(',k_string{i},'+',k_string{j},')*',k_string{i},'*',k_string{j},'.*',amp_string{i},'*',amp_string{j},'.*exp(1i*(',k_string{i},'+',k_string{j},')*x) '];
end
end
end
end
end
for j=1:length(w)
for i=1:length(w)
%%%% mixing_matrix(j,i) for j<i is dfg ie top right half
%%%% mixing_matrix(j,i) for j<i is sfg ie bottom left half
%%%% mixing_matrix(j,i) for j=i is SHG of ji = jj = ii is NW-SE
%%%% diagonal
if ismembertol(mixing_matrix(j,i),w,0.001)==1
if j<i
[check,place] = ismembertol(mixing_matrix(j,i),w,0.001);
dA_string_no_coeff{place,i} = ['+ ',k_string{i},'*',k_string{j},'*',amp_string{i},'*conj(',amp_string{j},')*exp(1i*(',k_string{i},'-',k_string{j},')*x) '];
end
if j>i
[check,place] = ismembertol(mixing_matrix(j,i),w,0.0001);
dA_string_no_coeff{place,i} = ['- ',k_string{i},'*',k_string{j},'*',amp_string{i},'*',amp_string{j},'.*exp(1i*(',k_string{i},'+',k_string{j},')*x) '];
end
if j==i
if ismembertol(mixing_matrix(j,i),w,0.0001)==1
[check,place] = ismembertol(mixing_matrix(j,i),w,0.001);
dA_string_no_coeff{i+j,i+j} = ['- 0.5*',k_string{i},'*',k_string{j},'.*',amp_string{i},'*',amp_string{j},'.*exp(1i*(',k_string{i},'+',k_string{j},')*x) '];
end
end
end
end
end
for i=1:length(dA_string)
dA_string_final{i} = ['(maxBeta/(2*',k_string{i},'))*(',dA_string{i,:},')*exp(-1i*',k_string{i},'*x);'];
end
for i=1:length(dA_string)
dA_string_no_coeff_final{i} = ['(maxBeta/(2))*(',dA_string_no_coeff{i,:},')*exp(-1i*',k_string{i},'*x);'];
end
fid = fopen([up_folder,'/CME_solver.m'],'w');
fprintf(fid, 'function [A] = CME_solver(maxBeta) \r\n\r\n');
for i=1:length(w)
k=acos(1-((w.^2).*param(3).*1e-12.*param(2).*1e-15)./(2.*(1-(w.^2).*param(3).*1e-12.*param(8).*1e-15)));
fprintf(fid, '%s = %0.4f\r\n',k_string{i},k(i));
end
fprintf(fid, '\r\n x=[0:%i] \r\n\r\n',N-1)
for i=1:length(dA_string_final)
if i==1
fprintf(fid, 'dA5 = @(x,A)[%s \r\n',dA_string_final{i});
elseif i==length(dA_string_final)
fprintf(fid, '%s] \r\n',dA_string_final{i});
else
fprintf(fid, '%s \r\n',dA_string_final{i});
end
end
fprintf(fid, '\r\n [x,A] = ode45(dA5, x, [')
for i =1:length(w)
if ismembertol(w(i),ws,0.001)==1
fprintf(fid, '%s ',As0)
elseif ismembertol(w(i),wp,0.001)==1
fprintf(fid, '%s ',Ap0)
else
fprintf(fid, '0 ')
end
end
fprintf(fid, '])')
fclose(fid)
fid = fopen([up_folder,'/CME_solver_no_coeff.m'],'w')
fprintf(fid, 'function [A] = CME_solver_no_coeff(maxBeta) \r\n\r\n')
for i=1:length(w)
k=acos(1-((w.^2).*param(3).*1e-12.*param(2).*1e-15)./(2.*(1-(w.^2).*param(3).*1e-12.*param(8).*1e-15)));
fprintf(fid, '%s = %0.4f\r\n',k_string{i},k(i))
end
fprintf(fid, '\r\n x=[0:%i] \r\n\r\n',N-1)
for i=1:length(dA_string_final)
if i==1
fprintf(fid, 'dA5 = @(x,A)[%s \r\n',dA_string_no_coeff_final{i})
elseif i==length(dA_string_final)
fprintf(fid, '%s] \r\n',dA_string_no_coeff_final{i})
else
fprintf(fid, '%s \r\n',dA_string_no_coeff_final{i})
end
end
fprintf(fid, '\r\n [x,A] = ode45(dA5, x, [')
for i =1:length(w)
if ismembertol(w(i),ws,0.001)==1
fprintf(fid, '%s ',As0)
elseif ismembertol(w(i),wp,0.001)==1
fprintf(fid, '%s ',Ap0)
else
fprintf(fid, '0 ')
end
end
fprintf(fid, '])')
fclose(fid)
%clear dA_string dA_string_final dA_string_no_coeff_final dA_string_no_coeff
Lg = param(3)*10^-12;
Ic = param(10)*10^-6;
Cg = param(2)*10^-15;
Cj = param(8)*10^-15;
phi0 = 2.067833848*10^-15;
Bl = 2*pi*Lg*Ic/phi0;
A=CME_solver(maxBeta);
V=(Lg*Ic/Bl).*w.*abs(A);
Z=sqrt(Lg/Cg);
I{harmonic} = V/Z;
A=CME_solver_no_coeff(maxBeta);
V=(Lg*Ic/Bl).*w.*abs(A);
Z=sqrt(Lg/Cg);
I_no_coeff{harmonic} = V/Z;
f_CME{harmonic}=w./(2*pi);
end
Thankyou for your consideration, it is much appreciated.
  2 Commenti
Sindar
Sindar il 13 Ott 2020
What is the actual error given?
Without looking too closely at your code, is it possible that the new version of the function has not been processed by Matlab, so it is using the old code where the size doesn't fit? This might be difficult to work around; is it possible you could generalize the function in such a way that it can take inputs to "update" rather than programmatically writing code?
Thomas Dixon
Thomas Dixon il 3 Dic 2020
Hi Sindar. Sorry for late response. I actually got this to work, but now I want to do exactly as you suggest. I dont want to have text files of pseudocode written and then read, as this seems fundementally stupid. (the worst kind of stupid). I want to somehow create the dA vector programaticaly based on the input A vector. I really am struggling with this at the moment. In the meantime I have created an app with default values which will calculate the mixing matrix and hopefully give insight into my problems.
You have to click 'Calculate' next to BetaL ( you may wish to click on the plot button as well to see my wonderful plotting)
Then go to the next tab ( I recommend that you choose MLine from the Cut-off determination tab and flick between All tones and Pump mediated for the Mixing mehcnaisms tab. The red light on Calculating tones means the script is running. I suggest for simplicity sticking with Pump mediated as it will be equal to or les than the number of entries in all tones.
Finally hit the Calculate Mixing Matrix buton and a table wil appear which has values corresponding with the values in the frequencies table (top right).
Whew.....
Now when you look at the mixing matrix table called mixingmatrix(i,j) then there is an obvious rule of correspondence between the numbers in Frequency=F(n) and MixingMatrix(i,j) the rule is:
for i=1:length(F)
for j=1:length(F)
if i>j
mixingmatrix(i,j)=F(i)-F(j)
end
if i<=j
mixingmatrix(i,j)=F(i)+F(j)
end
end
end
These mixing conditions and how theye are created define the relationship between A(n) and dA(n) which is always the same length of F(n). So with this rule and indexing I can work out the functional form of dA with regard to A. The problem I have is programtically giving it to matlab. I just don't know how to do it.
Thanks again for any help.

Accedi per commentare.

Risposte (0)

Categorie

Scopri di più su Startup and Shutdown in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by