Save each loop values, without overwriting

2016 4 27
1 个回答
更新时间：2016 4 27
21 次查看（30 天）

0 个投票

Hi. I think that my program is not saving each loop value for the Matrix Qabal which comes from Q, and for QQ neither, but I just want to solve the Q or Qabal, which it should be a matrix (22x24)
MND is a matrix (14x25) and MTB(which is inside the function has 22 rows).
for j=2:25
QQ=MND(:,j); %m3/s Vector containing fixed nodal discharges with each PF
[Y,Df,Q,VK1,Vn1]=funF2(H,QQ,MTB,VK,Vn,Vii,g,vi);
er=sum(abs(Y)); %Initial guess solution error
nit=0; %Iteration counter initiation
%Iterative solution search
while er>=eps && nit<=nma
VK=VK1;
Vn=Vn1;
H(Vii,j)=H(Vii)-Df\Y; %New extended H vector (solution vector)
[Y,Df,Q,VK1,Vn1,Qij]=funF2(H,QQ,MTB,VK,Vn,Vii,g,vi);
Qabal(Hores,j)=Q(Hores);
er=sum(abs(Y)); %New approximate solution error
nit=nit+1; %Iteration counter update
%Hores=Hores+1;
HH=H(:,2:end); %Only found H columns
end
end
I don't know If i have explained it well, maybe is needed the other part of the program to find it out, but I think that with that should be enough.
Thank you.

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

回答（1 个）

dpb
dpb 2016-4-27
编辑：dpb 2016-4-27

0 个投票

for j=2:25
QQ=MND(:,j-1); %m3/s Vector containing fixed nodal discharges with each PF
[Y,Df,Q,VK1,Vn1]=funF2(H,QQ,MTB,VK,Vn,Vii,g,vi);
The above will overwrite Q each pass with the results for that output argument from funF2 given the input using the j column of MND and the other inputs. We don't know what the dimensions are, presuming it is also a vector and you wish to save the various columns into an array Q instead, you would need
[Y,Df,Q(:,j),VK1,Vn1]=funF2(H,QQ,MTB,VK,Vn,Vii,g,vi);
instead and if that is so, you should "preallocate" Q as
Q=zeros(size(MTB,1),size(MND,1)-1); % preallocate
Similar holds for the other output array excepting you've written
Qabal(Hores,j)=Q(Hores);
which will load individual elements but you commented out the defining line that incremented the row index Hores. But, you can only have a single element to assign to a specific array location and it appears your function returns a vector so you'll need to (it would seem) also use
Qabal(:,j)=Q(Hores);
ERRATUM
Dang, missed an edit when pasted the line...intended to write
Qabal(:,j)=Q;
Of course, this raises the question of just how the iteration is supposed to be working--looks like maybe it's trying to update individual locations somehow, not the full vector but we don't know the intent just from code (particularly if the code doesn't work as intended :) ).

2 个评论
显示 无 隐藏 无

Jordi Escola
Jordi Escola 2016-4-27
I tried this before already but there is no way i can solve it. We are a group trying it but no solution...
Here is the "whole program", if you have time to check it, but don't worry if not, its quite long...
%node/length/diameter/e
MTB=[13,1,1000,0.380,0.0001,0,0;
1,2,1200,0.305,0.0001,0,0;
2,3,800,0.305,0.0001,0,0;
3,4,1200,0.255,0.0001,0,0;
2,5,1000,0.255,0.0001,0,0;
1,6,1200,0.255,0.0001,0,0;
1,7,1200,0.305,0.0001,0,0;
4,5,1000,0.205,0.0001,0,0;
5,6,1100,0.205,0.0001,0,0;
4,8,1200,0.255,0.0001,0,0;
5,9,1100,0.205,0.0001,0,0;
6,9,2000,0.205,0.0001,0,0;
6,10,1200,0.255,0.0001,0,0;
6,11,2000,0.255,0.0001,0,0;
7,11,2000,0.255,0.0001,0,0;
11,10,800,0.150,0.0001,0,0;
9,10,1000,0.205,0.0001,0,0;
8,9,800,0.150,0.0001,0,0;
8,12,1100,0.205,0.0001,0,0;
9,12,1000,0.205,0.0001,0,0;
10,12,2200,0.205,0.0001,0,0;
12,14,800,0.305,0.0001,0,0];
%height and Qout
MND=[Matrix 14x25, where i wanna take just columns 2 to 25]
Vii=[1,2,3,4,5,6,7,8,9,10,11,12]'; %Vector of nodes with unknown Hi
g=9.81; %N/kg Acceleration of gravity
vi=1.1e-6; %m2/s Fluid kinematic viscosity
Hores=[Vector of 24];
%Definition of computational parameters
eps=1e-6; %Maximum allowed error in the Fi
nma=30; %Maximum allowed iterative steps
%Definition of intermediate variables
nt=length(MTB(:,1)); %Number of pipes in the network
nn=length(MND(:,1)); %Total number of nodes in the network
ni=length(Vii); %Number of nodes with unknown Hi
npf=length(PF);
VK=zeros(nt,npf); %Vector containing the K value of each pipe
Vn=zeros(nt,npf); %Vector containing the n value of each pipe
%NEWTONS ITERATIVE METHOD
%Definition of the starting point
H=MND(:,1); %m Initial guess vector for H-eqs method (Hi=zi)
%Initial approximation to K and n of each pipe
f=0.02; %Initial guess for the friction factor
n=1.95; %Initial guess for exponent n
Q=zeros(size(MTB,1),size(MND,1)-1); % preallocate
for l=1:nt
if MTB(l,7)~=0 %Pipe element with fixed K and n
VK(l)=MTB(l,6);
Vn(l)=MTB(l,7);
else
VK(l)=8*f*MTB(l,3)/g/pi^2/MTB(l,4)^5; %Initial guess for K (eq. 19)
Vn(l)=n; %Initial guess for n
end
end
for j=2:25
QQ=MND(:,j-1); %m3/s Vector containing fixed nodal discharges with
[Y,Df,Q,VK1,Vn1]=funF2(H,QQ,MTB,VK,Vn,Vii,g,vi);
er=sum(abs(Y)); %Initial guess solution error
nit=0; %Iteration counter initiation
%Iterative solution search
while er>=eps && nit<=nma
VK=VK1;
Vn=Vn1;
H(Vii,j)=H(Vii)-Df\Y; %New extended H vector (solution
[Y,Df,Q(:,j),VK1,Vn1,Qij]=funF2(H,QQ,MTB,VK,Vn,Vii,g,vi);
Qabal(:,j)=Q(Hores);
er=sum(abs(Y)); %New approximate solution error
nit=nit+1; %Iteration counter update
%Hores=Hores+1;
HH=H(:,2:end); %Only found H columns
end
end
And the Function funF2:
%Auxiliar variables
nn=length(H); %Number of nodes in the network
nt=length(MTB(:,1)); %Number of pipes in the network
Y=zeros(nn,24); %Column vector containing all Fi functions
Df=zeros(nn); %Extended jacobian matrix initilization
Q=zeros(nt,24); %Column vector to contain the flow of each pipe
VK1=zeros(nt,24); %Column vector to contain new K values
Vn1=zeros(nt,24); %Column vector to contain new n values
%Iterative proces over all pipes
for l=1:nt
i=MTB(l,1);
j=MTB(l,2);
L=MTB(l,3); %m Pipe length
D=MTB(l,4); %m Pipe inner diameter
K=VK(l);
n=Vn(l);
Qij=sign(H(i)-H(j))*abs((H(i)-H(j))/K)^(1/n); %m3/s Flow from node i to node j
Q(l)=Qij;
%Computing the extended Fi terms
Y(i)=Y(i)+Qij;
Y(j)=Y(j)-Qij;
%Computing the extended jacobian matrix terms
Reij=4*abs(Qij)/pi/D/vi; %The actual Reynolds number
if Reij<=2300
dQdH=-1/K; %Laminar flow (eq.28)
else
dQdH=-abs((H(i)-H(j))/K)^(1/n-1)/n/K; %Turbulent flow (eq.26)
end
Df(i,j)=dQdH;
Df(j,i)=dQdH;
Df(i,i)=Df(i,i)-dQdH;
Df(j,j)=Df(j,j)-dQdH;
%Update of K and n for use in the next iteration
if MTB(l,7)~=0 %Pipe element with fixed K and n
VK1(l)=MTB(l,6);
Vn1(l)=MTB(l,7);
elseif Reij<=2300 %LAMINAR FLOW
VK1(l)=128*vi*L/pi/g/D^4; %Equation 17
Vn1(l)=1; %Equation 18
else %TURBULENT FLOW
Re=[0.95,1.05]*Reij; %The [Re1, Re2] points
f=1.325./log(MTB(l,5)/3.71/D+5.74./Re.^.9).^2; %Corresponding [f1, f2] points
b=log(f(1)/f(2))/log(Re(2)/Re(1)); %Equation 11
Q1=pi*D*vi*Re(1)/4;
a=f(1)*Q1^b; %Equation 12
VK1(l)=8*a*L/g/pi^2/D^5; %Equation 13
Vn1(l)=2-b; %Equation 14
end
end
Y=Y(Vii)+QQ(Vii); %m3/s Vector of Fi values for the ni nodes with unknown Hi
Df=Df(Vii,Vii); %The jacobian matrix
Thanks anyway, for the time that it took.
dpb
dpb 2016-4-27
Not sure what "this" refers to in the above but see the ERRATUM and comment following it to my previous Answer. Maybe the question raised there will lead to some enlightenment to resolve the issues.
As stated there, don't know any of us here can help without a more definitive statement of the problem trying to be solved, not just the code and a syntax issue with it.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Sparse Matrices 的更多信息

产品

标签

评论：

dpb
dpb
2016-4-27

Community Treasure Hunt

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

Start Hunting!

Translated by