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

3 views (last 30 days)
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.
  3 Comments
Ilham Hardy
Ilham Hardy on 3 Jul 2014
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
Anna on 4 Jul 2014
Edited: Cedric on 4 Jul 2014
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);

Sign in to comment.

Answers (0)

Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!