MATLAB Answers

Vector output of a function inside a for loop

8 views (last 30 days)
Meva
Meva on 25 Dec 2016
Edited: Meva on 26 Dec 2016
Hello,
So sorry for such a long code. I have a loop and a nested function whose input and output variables use the index of this loop. This index is n in the code.
When I debug, for n=1 u and p are in 1*101 and u are all 1 and p are all 0.
For n=2, u and p are 2*101, but the first rows are zero in u and p.
For n=3, u and p are 3*101, but the first and second rows are all zero in u and p.
So, Matlab forgets about the previous index. ii=101 in the code which is the second dimension of u and p, where first dimension is 3.
See attached please.
global ii
nn=3.;
ii=101;
b=1.;
dx=1./(ii-1.);
dt=0.001;
Nt=20001;
rmass=1.; rmomi=2.;
t=0.;
i_vect=(1:ii); x_vect=(i_vect-1)*dx;
dc=0.001; dp=0.001;
dckept=dc; dpkept=dp;
% Preallocation-----------------------------------------------------------%
...
%-------------------------------------------------------------------------%
%Initial conditions for H, A, HD, AD, cb, Anm1, ADnm1---------------------%
for n=1:nnsay
H(n)=1./nn; HD(n)=0.;
A(n)=0.;
abcd=(n-11.)^2.;
AD(n)=0.05*exp(-abcd);
A(nn)=0.; AD(nn)=0.;
cb(n)=1./nn;
Anm1=0.; ADnm1=0.;
if(n~=1)
Anm1=A(n-1);
ADnm1=AD(n-1);
end
end
%-------------------------------------------------------------------------%
for nt=1:Nt
% At the previous time step:
for n=1:nn
for i=1:ii
x=(i-1.)*dx;
ub(n,i)=(cb(n)-HD(n)*(x-b/2.)-(ADnm1-AD(n))/2.*(x-b/2.)^2)...
/(H(n)+(Anm1-A(n))*(x-b/2.));
pb(n,i)=0.;
end
end
for n=1:nn
Hbar(n)=H(n);
HDbar(n)=HD(n);
Abar(n)=A(n);
ADbar(n)=AD(n);
c(n)=cb(n);
end
flag3=0;
A(nn)=0.;
AD(nn)=0.;
for lmn=1:8
mmm=1;
pend=pb(1,ii); % At the first time step, pend=0;
flag2=10;
while (flag2==10)
m=1;
while (flag3==0)
for n=1:nn
kkk=1;
k = 1;
nm=n;
[this_u,this_p]=uqp(ii,nm, dx,dt,ub,c,HD,b,AD,H,A);
u=[u; this_u];
p=[p; this_p];
st(k)=p(n,ii)-pend;
k=k+1;
if(k==2)
c(n)=c(n)+dc;
end
if(k==3)
den=st(1)-st(2);
c(n)=(c(n)*st(1)-(c(n)-dc)*st(2))/den;
end
if (k==4)
kkk=kkk+1;
if (kkk==2)
flag3=1;
end
end
dc=dckept;
end
end
ss=0.;
for n=1:nn
ss=ss+c(n);
end
rt(m)=ss-1.;
m=m+1;
if (m==2)
pend=pend+dp;
end
if (m==3)
deno=rt(1)-rt(2);
pend=(pend*rt(1)-(pend-dp)*rt(2))/deno;
end
if (m==4)
mmm=mmm+1;
if (mmm==2)
m=1;
else
flag2=11;
end
end
dp=dpkept;
end
end
for n=1:nn
cb(n)=c(n);
for i=1:ii
ub(n,i)=u(n,i);
pb(n,i)=p(n,i);
end
end
t = t + dt;
end
which is the main code. The nested function is:
function [u,p]=uqp(ii,nm,dx,dt,ub,c,HD,b,AD,H,A)
n=nm;
Anm1=0.;
ADnm1=0.;
ub(n,:)
if(n~=1)
Anm1=A(n-1);
ADnm1=AD(n-1);
end
for i=1:ii
x=(i-1.)*dx;
u(n,i)=(c(n)-HD(n)*(x-b/2.)-(ADnm1-AD(n))/2.*(x-b/2.)^2)/(H(n)+(Anm1-A(n))*(x-b/2.));
end
p(n,1)=0.5*(1.-ub(n,1)*u(n,1));
q(n,1)=0.;
for i=2:ii
q(n,i)=q(n,i-1)-dx*(u(n,i-1)-ub(n,i-1))/dt;
p(n,i)=0.5*(1.-ub(n,i)*u(n,i))+q(n,i);
end
return
end
Any idea will highly be appreciated! In attached, main.m is main code which uses secant method double times. uqp.m and HA.m are nested functions.

  2 Comments

John BG
John BG on 25 Dec 2016
why do i have the impression that such long string of input variables, perhaps one fell off and by not passing the variable MATLAB has to assume or ignore?
would you please be so kind to supply a chunk of working code with variables. It hasn't got to be actual code, or the actual values in variables, but it's highly appreciated by readers to have some kind of code with variables, you only supplied part of the algorithm.
When asking, a bit too much of information is usually better than a few too less, don't you agree?
Meva
Meva on 25 Dec 2016
Yes, of course I do agree! I just thought that it would be messy if I pasted the whole code. I will edit my question accordingly then! Thanks @Jofn BG

Sign in to comment.

Answers (1)

Image Analyst
Image Analyst on 25 Dec 2016
Edited: Image Analyst on 25 Dec 2016
You need to pass p in if you want it to retain all values, or else have the function return just the newest row, and use a temporary p called "this_p" and append it to the "master" p like this:
[u, this_p] = uqp(ii, nm, dx, dt, ub, c, HD, b, AD, H, A);
p = [p; this_p]; % Append new row to the bottom of the existing p.

  4 Comments

Show 1 older comment
Meva
Meva on 25 Dec 2016
If I work with this_u and this_p, it increases 3 times in every loop I think.
Image Analyst
Image Analyst on 26 Dec 2016
I copied and pasted, added a function line to the main program and it did not run. Please attach the whole m-file with the paper clip icon.
Meva
Meva on 26 Dec 2016
Thanks so much @Image Analyst . Yes, I cropped my code to make it less messy.
Now, I attach my all codes; main.m is main code which uses secant method double times.
uqp.m and HA.m are nested functions to calculate u, q, p and H and A respectively.
When I run my main program, u does not change over x but it should have changed. Another issue is that I think H(n), HD(n), AD(n), A(n) did not go inside properly as an input of uqp.m function.
Thanks so much for your help in advance.

Sign in to comment.

Sign in to answer this question.


Translated by