Convergence program stuck at a point.

Hi all,

Below is the code when I am trying to get convergence of a function by varying a variable which dictates it. The sum of variables is a fixed constant. i.e a+b = a_initial. The issue I am facing is after some iterations a_old=a_older and b_old=b_older and the new values don't update. I know that there is a fix for this but somehow I am unable to fix it. Any help would be appreciated.

 a_initial = M_f_b_i;
a_older = M_f_b_i; % initial guess for a
b_older = 0;% intial guess for b
tolerance = 1e-1; % tolerance for the convergence
b_old = (a_older+b_older)/2;
a_old = M_f_b_i-b_old;
kb = get_curvature_bending(a_old,Eb);
kr = get_curvature_rotating(b_old,Er);
i=1;
while abs(kb-kr)>tolerance
    if kb>kr % 
        b_new = (b_old+a_old)/2;
        a_new = (a_initial-b_new);
        kb = get_curvature_bending(a_new,Eb);
        kr = get_curvature_rotating(b_new,Er);
        i=i+1;
        a_older = a_old;
        b_older = b_old;
        a_old = a_new;
        b_old = b_new;
    else    
        b_new = (b_old+b_older)/2;
        a_new = (a_initial-b_new);
        kb = get_curvature_bending(a_new,Eb);
        kr = get_curvature_rotating(b_new,Er);
        i=i+1;

        a_older = a_old;
        b_older = b_old;
        a_old = a_new;
        b_old = b_new;
    end
    sprintf('Iteration number: %d',i)
end
end_curvature = kb;
Moment_f_b = a_new;
Moment_f_r = b_new;
end

Answer 1

[Roger Stafford]

I think you are confusing things with your names of ‘old’ and ‘older’. What you need are the concepts of “too high” and “too low”. Since you are restricting the sum of ‘a’ and ‘b’ to be a given constant you might as well do things entirely in terms of, say, ‘a’, and alway use ‘b’ as being equal to that constant minus ‘a’. Next you have to decide in advance which way kb-kr will be when a is too high and which way it will be when a is too low. If it will be positive when a is too high and negative when a is too low, you can do something like this:

   set a_too_high so that it is certainly too high
   set a_too_low so that it is certainly too low
   set kb-kr beyond tolerance
   while abs(kb-kr)>tolerance
     a_test = (a_too_high+a_too_low)/2
     get b_test and compute kb-kr
     if kb-kr>0  % a_test is too high, move a_too_high down
       a_too_high = a_test;
     else        % a_test is too low, move a_too_low up
       a_too_low = a_test
     end
   end

This way you should eventually find some value of a in between the original too high and too low values when it is just right.