Shorten computation time for solving recursive equation

Greetings.

I'm trying to solve a recursive symbolic equation. a is the variable I'm solving for and c, d and e are constants. It goes like this:

var = c+d*a/(1-e*a);

for 1:N

var=c+d*(var)/(1-e*var);

end

Then I set var == 0 and numerically solve for a. It takes a very long time for Matlab to solve when N increases, even with a decent initial guess. Is there any way to shorten the computation time?

Thanks.

[Walter Roberson]

What sort of order of magnitude is N ?

And how are you doing the numeric solving?

The small test I just did with N = 10 generated a var with powers of a above 1300 (display got too wide at that point); vpasolve() of that without any initial condition would try to find all 1300+ roots.

What sort of magnitude are c and d? I just tested with 2 and 9, and the terms involved coefficients up to about 10^1382 -- far too large for double precision to be useable.

[Jerry Yeung]

Thanks for the reply. Apologies but I've made a mistake in my initial post. Please see the amended version. c, d and e are all complex with their magnitudes smaller than 1. N is on the order of 10e3.

Answer 1

[Walter Roberson]

N = 100;
c = complex(randn,randn); d = complex(randn,randn); e = complex(randn,randn);
c = c./abs(c).*rand; d = d./abs(d).*rand; e = e./abs(e).*rand;
disp([c,d,e])

  -0.5490 - 0.1728i  -0.6121 - 0.4458i   0.4991 - 0.5406i

tic
syms a

var = c+d*a/(1-e*a);
for K = 1:N
    var = simplify(c+d*(var)/(1-e*var));
end
toc

Elapsed time is 29.084258 seconds.

tic
A = solve(var); 
toc

Elapsed time is 0.356109 seconds.

[Jerry Yeung]

In your example the code still took 30 seconds to run. I have a range of c,d and e constants to sweep through and this can easily take hours. Can you suggest ways to shorten the runtime?

[Walter Roberson]

There is a possibility that working purely numerically might work. That is, instead of calculating the formula and solve() it, instead accept a numeric trial a, run through the loop numerically, and fsolve() the whole thing or something like that.

That said, I seem to be having difficulty getting the numeric form to work, whereas the symbolic form works fine (see the cross-check.)

This has a few potential explanations:

• there might not be enough precision in doubles to resolve the equations; you might need to balance expressions to more than 16 digits
• the simplify() step that reduces the formula to make the solve() fast, might be implicitly gaining precision by reordering operations -- fewer terms evaluated with a specific a might drastically reduce the round-off problems
• there might be sub-expressions underflowing to 0 when done numerically. If you look at the function being solved for 0 in the symbolic version, you will find it involves terms with values on the order of 1E+3600. We are not seeing inf in the numeric versions of the calculation, but perhaps some of the terms underflow below 1e-308 and we lose the balance.

outer

  -0.5495 - 0.5652i  -0.2913 + 0.9064i   0.5610 - 0.2352i


Initial point is a local minimum that satisfies the constraints.

Optimization completed because at the initial point, the objective function is non-decreasing 
in feasible directions to within the value of the optimality tolerance, and 
constraints are satisfied to within the value of the constraint tolerance.

numeric determined a was 2 in 0.29994 seconds
Residue:     0.6094

symbolic determined a was 2.063262916322103 in 29.26626 seconds

crosscheck_sym_num = 0.6096

crosscheck_sym_sym = 

0

function outer
    rng(655321)
    
    N = 100;
    c = complex(randn,randn); d = complex(randn,randn); e = complex(randn,randn);
    c = c./abs(c).*rand; d = d./abs(d).*rand; e = e./abs(e).*rand;
    disp([c,d,e])

    tic
    guess = 2;
    [Anum, fval] = fmincon(@calculate_a_numerically, guess);
    Anum_time = toc;
    fprintf('numeric determined a was %.16g in %.5f seconds\nResidue: ', Anum, Anum_time);
    disp(fval)

    tic
    Asym = calculate_a_symbolically();
    Asym_time = toc;
    fprintf('symbolic determined a was %.16g in %.5f seconds\n', Asym, Asym_time);

    crosscheck_sym_num = calculate_a_numerically(double(Asym))
    crosscheck_sym_sym = calculate_a_symbolically(Asym)

    function A = calculate_a_symbolically(a)    %nested function!
        if nargin < 1;
            a = sym('a');
            dosolve = true;
        else
            dosolve = false;
        end
        
        var = c+d*a/(1-e*a);
        for K = 1:N
            var = simplify(c+d*(var)/(1-e*var));
        end
        if dosolve
            A = solve(var);
        else
            A = var;
        end
    end
    
    function A = calculate_a_numerically(a)
        var = c+d*a/(1-e*a);
        for K = 1:N
            var = c+d*(var)/(1-e*var);
        end
        %A = [real(var), imag(var)];
        A = norm(var);
    end
end

[Jerry Yeung]

Apologies as it took me a while to verify this. Indeed, the numerical approach doesn't always return a solution, despite the symbolic approach yielding a perfect valid solution even for small N. I suppose I'll have to code my own iterative approach for large N then. What Matlab functions would help with this?