Vector output of a function inside a for loop

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.

[John BG]

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]

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

Answer 1

[Image Analyst]

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.

[Meva]

Thanks @Image Analysts. It could not seem to work to me I should say.I will paste my code with more info.

[Meva]

If I work with this_u and this_p, it increases 3 times in every loop I think.

[Image Analyst]

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]

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.