How can I write the following Algorithm in MATLAB? What is wrong with my code?

Hello,

We have an algorithm like:

Step 1: Set ε, ρ_0,t_1. Set i=1

Step2: If a≤ε, stop.

Step 3: Calculate ∆_i. Set t_{i}^{+}=t_i+∆_i

Step 4: Calculate ρ_i.

Step 5: If b=0, go to Step 6. Otherwise, set t_{i+1}=t_{i}^{+} and go to Step 7.

Step 6: If ρ_i ≥ρ_0, set t_{i+1}=t_{i}^{+}, do something go to Step 8. Otherwise do another thing and go to Step 3 (inner cycle).

Step 7: If ρ_i ≥ρ_0 do something and go to Step 8. Otherwise do another thing and go to Step 8.

Step 8: Set i=i+1 and go to Step 2. (outer cycle)

I have written a code like this:

Set ε, ρ_0,t_1, maxiter. Calculate a, Set i=1

While (a> ε & i<maxiter)

          Calculate ∆_i
          Set t_{i}^{+}=t_i+∆_i
          Calculate ρ_i.
          Calculate b.
If b=0 
   If ρ_i ≥ρ_0
         t_{i+1}=t_{i}^{+}
         i=i+1
else
        while (ρ_i <ρ_0  & b=0)
               Calculate ∆_i
               Set t_{i}^{+}=t_i+∆_i
               Calculate ρ_i.
               Calculate b.
        end   (out of inner cycle)
If ρ_i ≥ρ_0  
            t_{i+1}=t_{i}^{+}
            i=i+1
            do something
else
                     t_{i+1}=t_{i}^{+}
              i=i+1
            do another thing
end
end
else  (beginning if)
  t_{i+1}=t_{i}^{+}
If ρ_i ≥ρ_0  
    do something
      i=i+1
else
       do another thing 
       i=i+1
end
calculate a
end (out of  outer cycle)

The code is working but it is not working as it should work, I could not find my mistake in the code. Do you have any suggestions for the code? Especially for the loop (inner cycle) between Step 3 and Step 6.

[Geoff Hayes]

Anna - what do you mean by the code is working but it is not working as it should work? Is the output unexpected given the input? Is an error thrown?

It might be more helpful if you could attach the code to your question so that it could be compared (more easily) to the algorithm. While the above pseudo-code is helpful, it is not formatted cleanly (and is incomplete) making it difficult to determine where it might deviate from the algo.

[Ilham Hardy]

To be working and not to be working.. :D

Could you be more specific to describe what you see and what you got?

It is also important to put error messages (if there any) to help us help you.

Regards, IH

(PS. no pun intended)

[Anna]

Here is my Matlab function to find zero of a nonlinear function f(z). It does not converge to solution the output of the function should go to zero (it should converge at least locally, when we begin in a neighborhood of exact solution) but it does not.

Could you find my mistake in the code? Do I apply Algorithm above in a wrong way?

 function dz=f1zero(z) 

 %input: z is a 5x1 vector

 %WE WANT TO SOLVE A NONLINEAR SYSTEM OF EQUATIONS: %H(z)=0, see function below f(z) below.

 %Exact solution is: (found by fsolve) %z=[-0.0953101798041584;0.0953101798044882;1.09999999999948;1.00000000000018;0.999999999996862];

 % Step 1 

 format long g 

 %set parameters 

 h=10; 
 rhozero=0.1; 
 gtheta=0.001;

 %for termination criteria 

 tol=1e-6; 

 itermax=20; %maximum number of iterations 

 %iteration count 
 k=1; 

 Hg=f(z); %nonlinear system of equations we want to solve H(z)=0 

 Wg=jacf(z); %Jacobian of the system H(z) 

 grad=Wg'*Hg; %gradient of F(z)=1/2||H(z)||^2

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BEGİNNİNG OF OUTER CYCLE

 while (norm(grad,2)>tol & k<itermax) 

 k; 
 Hg=f(z); 
 Wg=jacf(z); 
 grad=Wg'*Hg; 
 VV=Wg'*Wg+((1/h)*eye(5)); 
 d=-grad\VV;

 if k==1 FF(1,:)=ctheta(Hg); end

 zplus=z+d'; 
 qkd=(0.5)*(norm(Hg+Wg*d',2)^2); 

 Fn=Thetaf(z); %see below

 Hg=f(zplus);                %system at zplus
 Fnz=Thetaf(zplus);          %norm at zplus
 rho=(Fn-Fnz)/(Fn-qkd);      %calculate rho_k 
 thzp=ctheta(Hg);

 [a,fi]=Accept(FF,thzp); FF=fi;

     if a==0                 
                if rho>=rhozero   
                   z=zplus;
                   h=2*h;
                   k=k+1;
                else   
                   it=0;
            while (rho<rhozero & a==0)   
               it=it+1;
               h=(0.5)*h;
               VV=Wg'*Wg+((1/h)*eye(5));
               d=-grad\VV;
               znplus=z+d'; 
               Hg=f(znplus);                
               thznp=ctheta(Hg);
               [a,fi]=Accept(FF,thznp);   
               Fnz=Thetaf(znplus);            
               rho=(Fn-Fnz)/(Fn-qkd) ;         
             end %end of while inner cycle         
          if rho>=rhozero
              z=znplus;
              h=2*h;
              k=k+1;
          else
              z=znplus;    
              h=0.5*h;
              k=k+1;
              FF=fi;
          end
         end %if rho>=rhozero   
     else    %a=1 
         z=zplus;  
        if rho>=rhozero
              h=2*h;
              k=k+1;
        else 
              [a,fi]=Accept(FF,thzp);
              FF=fi;   
              h=(0.5)*h;
              k=k+1;
          end
     end       %end of beginning if  

 Hg=f(z);
 Wg=jacf(z);
 grad=Wg'*Hg;
 NORM=norm(grad,2)

 end %end of while outer cycle

 %Outputs: func. values at final iterate

 Hg=f(z); Wg=jacf(z); grad=Wg'*Hg; dz=norm(grad,2)

 %------------------------------------ 

 function d=f(z) 

 %Separating z=(x,u,y,w) 

 n=2;m=1;p=1;q=1;

 x=z(1:n,:); u=z(n+1:n+p,:); y=z(n+p+1:n+p+m*p,:); w=z(n+p+m*p+1:n+p+m*p+q*p,:);

 d=[1.21*exp(x(1))-u*exp(x(1)+x(2)); exp(x(2))-u*exp(x(1)+x(2)); sqrt(u^2+(y-exp(x(1)+x(2)))^2)-u+y-exp(x(1)+x(2)); 1-2*w*y+w; sqrt(w^2+(y^2-y)^2)-w+y^2-y ];

 %-----------------------------------------------

 function jd=jacf(z)

 %finds Jacobian of f

 n=2;m=1;p=1;q=1;

 x=z(1:n,:); u=z(n+1:n+p,:); y=z(n+p+1:n+p+m*p,:); w=z(n+p+m*p+1:n+p+m*p+q*p,:);

 tol=1e-6; if abs(u)<tol & abs(y-exp(x(1)+x(2)))<tol lam=1; gam=-1; else lam=(y-exp(x(1)+x(2)))/sqrt(u^2+(y-exp(x(1)+x(2)))^2)+1; gam=u/(sqrt(u^2+(y-exp(x(1)+x(2)))^2))-1; end

 if abs(w)<tol & abs(y^2-y)<tol alp=1; bet=-1; else alp=(y^2-y)/sqrt(w^2+(y^2-y)^2)+1; bet=(w/sqrt(w^2+(y^2-y)^2))-1; end jd=[1.21*exp(x(1))-u*exp(x(1)+x(2)) -u*exp(x(1)+x(2)) -exp(x(1)+x(2)) 0 0; -u*exp(x(1)+x(2)) exp(x(2))-u*exp(x(1)+x(2)) -exp(x(1)+x(2)) 0 0; -lam*exp(x(1)+x(2)) -lam*exp(x(1)+x(2)) gam lam 0; 0 0 0 -2*w 1-2*y; 0 0 0 alp*(2*y-1) bet];

 %--------------------------------------------------- 

 function T=Thetaf(z) 

 %evaluates theta(z)=(1/2)*H(z)'*H(z) 

 T=.5*f(z)'*f(z);