Shorten computation time for solving recursive equation

Question

Jerry Yeung 2021-3-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/775462-shorten-computation-time-for-solving-recursive-equation

编辑： Walter Roberson 2021-3-24

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.

2 个评论
显示无隐藏无

Walter Roberson 2021-3-17

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 2021-3-17

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 2021-3-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/775462-shorten-computation-time-for-solving-recursive-equation#answer_650422

编辑：Walter Roberson 2021-3-17

在 MATLAB Online 中打开

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Walter Roberson 2021-3-18

编辑：Walter Roberson 2021-3-24

在 MATLAB Online 中打开

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 2021-3-24

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?

请先登录，再进行评论。

Shorten computation time for solving recursive equation

2 个评论
显示无隐藏无

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Shorten computation time for solving recursive equation

2 个评论 显示 无隐藏 无

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

3 个评论
显示 1更早的评论隐藏 1更早的评论