Iterating formula until value converges

So the problem is i would like to use a while loop to iterate an equation by substituting values calculated in the last iteration into the next iteration until the value epsilon_conv converges to some small number such as 0.001. Below is the code so far, i am able to find the first iteration ok but then i'm not sure how to sub the new values in and use the equation for epsilon_conv with the old and new beta values to find convergence. Thanks in advance.

clear all;
syms Sy Pa b h Pb l
f=(Sy-(Pa/(b*h)))-(6*Pb*l/(b*h^2));  % master equation

dSy=diff(f,Sy);          % f diffirentiated wrt Sy
dPa=diff(f,Pa);          % f diffirentiated wrt Pa
db=diff(f,b);            % f diffirentiated wrt b
dh=diff(f,h);            % f diffirentiated wrt h
dPb=diff(f,Pb);          % f diffirentiated wrt Pb
dl=diff(f,l);            % f diffirentiated wrt l

Sy=4000;               % values used in first iteration
Pa=2500;
b=2;
h=5;
Pb=900;
l=16;

Sys=400;               % values used in first iteration
Pas=250;
bs=0.2;
hs=0.5;
Pbs=90;
ls=1.6;


epsilon=1;
epsilon_conv=0.001;
i=1;
while epsilon>epsilon_conv
    
    if i=1
        oldbeta=1;
    else
        
        
        
        DSy=subs(dSy);           % Substituting values into diffirentiated equations
        DPa=subs(dPa);
        Db=subs(db);
        Dh=subs(dh);
        DPb=subs(dPb);
        Dl=subs(dl);
        fnorm=subs(f);
        
        beta1=fnorm/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);   % defining beta
        
        a11=-(DSy*Sys)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
        a21=-(DPa*Pas)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
        a31=-(Db*bs)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
        a41=-(Dh*hs)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
        a51=-(DPb*Pbs)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
        a61=-(Dl*ls)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
        
        Sy=Sy+(beta1*Sys*a11);   %%% new value of Sy etc. to be used in calculation of next iteration
        Pa=Pa+(beta1*Pas*a21);
        b=b+(beta1*bs*a31);
        h=h+(beta1*hs*a41);
        Pb=Pb+(beta1*Pbs*a51);
        l=l+(beta1*ls*a61);
        
        Sys=(x12-Sy)/Sys;   %% new value of Sys etc. to be used in calculation of next iteration
        Pas=(x22-Pa)/Pas;
        bs=(x32-b)/bs;
        hs=(x42-h)/hs;
        Pbs=(x52-Pb)/Pbs;
        ls=(x62-l)/ls;
        
        epsilon_conv=(oldbeta-newbeta)/newbeta  %% number that must converge to 0.001

Answer 1

[Ameer Hamza]

Overall the idea of your code is correct, but there are a few syntax issues.

The line

if i=1

is incorrect, use i==1 as condition.

I guess you just want to run the line old_beta = 1 once in first iteration. If yes, then you should change the variable i so that its is not repeated again

if i==1
    oldbeta=1;
    i = 2;
else

Also in these lines

Sys=(x12-Sy)/Sys;   %% new value of Sys etc. to be used in calculation of next iteration
Pas=(x22-Pa)/Pas;
bs=(x32-b)/bs;
hs=(x42-h)/hs;
Pbs=(x52-Pb)/Pbs;
ls=(x62-l)/ls;

valuesof x12, x22, x32, ..., x62 are not defined. Define them before use.

Since you are using symbolic variables, your code will run slow. If fast processing is required, study about matlabFunction() and function handless.

[ab123]

The code running slowly shouldn't be a problem really. The undefined x values were an accident as id changed pard of the code and not changed that bit. This is more what ive got now but it gives an error on line 49 saying Dsy is undefined, and i'm not really sure how to progress

clear all;
syms Sy Pa b h Pb l
f=(Sy-(Pa/(b*h)))-(6*Pb*l/(b*h^2));  % master equation

dSy=diff(f,Sy);          % f diffirentiated wrt Sy
dPa=diff(f,Pa);          % f diffirentiated wrt Pa
db=diff(f,b);            % f diffirentiated wrt b
dh=diff(f,h);            % f diffirentiated wrt h
dPb=diff(f,Pb);          % f diffirentiated wrt Pb
dl=diff(f,l);            % f diffirentiated wrt l

Sy=4000;               % values used in first iteration
Pa=2500;
b=2;
h=5;
Pb=900;
l=16;

Sys=400;               % values used in first iteration
Pas=250;
bs=0.2;
hs=0.5;
Pbs=90;
ls=1.6;


epsilon=1;
epsilon_conv=0.001;
i=1;
while epsilon>epsilon_conv
    
    if i==1
        oldbeta=1;
        
    else
        
        DSy=eval(subs(dSy));           % Substituting values into diffirentiated equations
        DPa=eval(subs(dPa));
        Db=eval(subs(db));
        Dh=eval(subs(dh));
        DPb=eval(subs(dPb));
        Dl=eval(subs(dl));
        fnorm=eval(subs(f));
        
        beta1=fnorm/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);   % defining beta
        i=i+1;
    end
    
    a11=-(DSy*Sys)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
    a21=-(DPa*Pas)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
    a31=-(Db*bs)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
    a41=-(Dh*hs)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
    a51=-(DPb*Pbs)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
    a61=-(Dl*ls)/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);
    
    x12=Sy+(beta1*Sys*a11);
    x22=Pa+(beta1*Pas*a21);
    x32=b+(beta1*bs*a31);
    x42=h+(beta1*hs*a41);
    x52=Pb+(beta1*Pbs*a51);
    x62=l+(beta1*ls*a61);
    
    u12=(x12-Sy)/Sys;
    u22=(x22-Pa)/Pas;
    u32=(x32-b)/bs;
    u42=(x42-h)/hs;
    u52=(x52-Pb)/Pbs;
    u62=(x62-l)/ls;
    
    Sy=x12;    % new value of Sy etc. to be used in calculation of next iteration
    Pa=x22;
    b=x32;
    h=x42;
    Pb=x52;
    l=x62;
    
    Sys=u12;  % new value of Sys etc. to be used in calculation of next iteration
    Pa=u22;
    bs=u32;
    hs=u42;
    Pb=u52;
    ls=u62;
    
    
    epsilon=(oldbeta-beta)/beta;  %% number that must converge to 0.001
end

[Ameer Hamza]

abreg123, in first iteration, the condition for if block become true and only this lines runs

oldbeta=1;

Following lines don't run in first iteration

DSy=eval(subs(dSy));           % Substituting values into diffirentiated equations
DPa=eval(subs(dPa));
Db=eval(subs(db));
Dh=eval(subs(dh));
DPb=eval(subs(dPb));
Dl=eval(subs(dl));
fnorm=eval(subs(f));

beta1=fnorm/sqrt((DSy*Sys)^2+(DPa*Pas)^2+(Db*bs)^2+(Dh*hs)^2+(DPb*Pbs)^2+(Dl*ls)^2);   % defining beta
i=i+1;

so the value of Dsy is not defined and it gives error on first iteration.

Also note that the value of the variable i will never change in your code because the initial value of i is 1, the if-condition becomes true. But the value of i does not change, because the else block never run. Also, note that beta is not defined in your code. Can you explain the program logic, maybe I can suggest an easier approach.