Evaluate a changing expression in loop

Hi! I am writing a code that generates a function f in a loop. This function f changes in every loop, for example from f = x + 2*x to f = 3*x^2 + 1 (randomly), and I want to evaluate f at different points in every loop. I have tried using subs, eval/feval, matlabFunction etc but it is still running slowly. How would you tackle a problem like this in the most efficient way?

Answer 1

[Stephen23]

It all depends on what you mean by "random": surely there must be some limit to the power or coefficient, and a limit to the number of addends? What distribution do these values have?

The two examples that you show both have two addends, and each addend is a low order power of x with a coefficient. If these are indicative of all of your functions, then why not use MATLAB effectively, instead of trying to fight it with evil eval? MATLAB works best with matrices, so lets try using them! Simply generate all of the coefficients, powers, etc into matrices, and calculate the function once, when required. This will be much faster and more efficient than generating "random" functions in a loop and evaluating them.

Here is an example showing twenty (5x4) independent, random functions, all calculated simultaneously by using 5x4 arrays for the function coefficients and powers.

>> cof = randi(4,5,4,2)-1; % all 20+20 coefficients
>> pwr = randi(4,5,4,2)-1; % all 20+20 powers
>> F = @(x)cof(:,:,1).*x.^pwr(:,:,1) + cof(:,:,2).*x.^pwr(:,:,2);
>> F(1) % all 20 functions at x=1
ans =
     2     1     6     0
     3     4     3     3
     3     5     1     3
     4     3     3     2
     1     0     1     0
>> F(2) % all 20 functions at x=2
ans =
     9     2    24     0
     6    16     6    24
     6    20     4    12
     7    16     6     8
     2     0     4     0

Lets have a look in more detail at the first element of the output. The coefficients and powers of its two addends are:

>> cof(1,1,:)
ans(:,:,1) =
     1
ans(:,:,2) =
     1
>> pwr(1,1,:)
ans(:,:,1) =
     0
ans(:,:,2) =
     3

Note that the pwr of zero makes this addend a constant, equal to cof(1,1,1), because x^0 is one. So for x==2, the first element is thus:

cof(1,1,1)*2^pwr(1,1,1) + cof(1,1,2)*2^pwr(1,1,2)

which is:

1*2^0 + 1*2^3

which is:

1*1 + 1*8

which is equal to

9

the same as the matrix shows for the first element (at x==2).

[Kristoffer Lindvall]

Thank you very much for your answer! I think that writing randomly was a mistake because the functions that I generate are not random. I wrote that as a thought experiment. I have to evaluate the function f at a point x0 before continuing to the next step. I will take into consideration what you wrote. I might be going about this in a complicated manner.

[Kristoffer Lindvall]

The code below is my solution and it is one loop

    x = sym('x',[2,1]);
    f = [x(1)-x(1)cos(x(2)), x(2)-3x(2)^2*cos(x(1))];
    J = jacobian(f,x);
    x0 = [2,1];
    N=length(x0); % Number of equations

    %%Transform into string
    fstr = map2mat(char(f));
    Jstr = map2mat(char(J));

    % replace every occurence of 'xi' with 'x(i)'
    Jstr = addPar(Jstr,N);
    fstr = addPar(fstr,N);

    x = x0;

    phi0 = eval(fstr)
    J = eval(Jstr)

    function str = addPar(str,N)

    % pstr = addPar(str,N)
    % Transforms every occurence of xi in str into x(i)
    % N is the maximum value of i
    % replace every occurence of xi with x(i)
    % note that we do this backwards to avoid x10 being
    % replaced with x(1)0

    for i=N:-1:1
    is = num2str(i);
    xis = ['x' is];
    xpis = ['x(' is ')'];
    str = strrep(str,xis,xpis);
    end

    function r = map2mat(r)

    % MAP2MAT Maple to MATLAB string conversion.
    % Lifted from the symbolic toolbox source code
    % MAP2MAT(r) converts the Maple string r containing
    % matrix, vector, or array to a valid MATLAB string.
    %
    % Examples: map2mat(matrix([[a,b], [c,d]]) returns
    % [a,b;c,d]
    % map2mat(array([[a,b], [c,d]]) returns
    % [a,b;c,d]
    % map2mat(vector([[a,b,c,d]]) returns
    % [a,b,c,d]
    % Deblank.

    r(findstr(r,' ')) = [];

    % Special case of the empty matrix or vector
    if strcmp(r,'vector([])') | strcmp(r,'matrix([])') | ...
    strcmp(r,'array([])')
    r = [];
    else
    % Remove matrix, vector, or array from the string.
    r = strrep(r,'matrix([[','['); r = strrep(r,'array([[','[');
    r = strrep(r,'vector([','['); r = strrep(r,'],[',';');
    r = strrep(r,']])',']'); r = strrep(r,'])',']');
    end
    end

[Stephen23]

What you are attempting is going to be a very slow way of achieving your aims. Mucking around concatenating strings and then relying on string evaluation is especially slow and leads to buggy code that is impossible to debug, as is explained here:

http://www.mathworks.com/matlabcentral/answers/196952-how-to-create-new-variables-in-batches-with-strcat

I would suggest that you learn about more efficient ways of programming in MATLAB, and taking advantage of its array handling abilities. I notice that all of your questions on this forum are basically "I write slow code using eval... help me": if you want faster, more efficient code, think about what MATLAB does best and fast: processing of numeric arrays. Learning to think in MATLAB terms means learning some new ways of thinking about algorithms, code and data:

http://www.mathworks.com/matlabcentral/answers/228557-experts-of-matlab-how-did-you-learn-any-advice-for-beginner-intermediate-users

Perhaps you could actually tell us what your aim is, what your are trying to achieve, and we can suggest efficient ways of solving your task.

[Kristoffer Lindvall]

Yes! You are absolutely correct, I will take that criticism to heart. I started using eval, begrudgingly I might add, only because it went a factor 6 times faster than subs, and twice as fast as matlabFunction.

My goal is to use Chebyshev polynomials to approximate the solution to ODEs and PDEs. I am using a semi implicit root solver that solves the system of equations in every time interval. The SIR algorithm however needs a symbolic expression of the function f and the jacobian J, which I need to evaluate at the beginning of every time interval. I have struggled to find a way to do this numerically, but it seems as if I am stuck with doing it analytically. The algorithm for SIR can be found at http://uu.diva-portal.org/smash/record.jsf?pid=diva2%3A217238&dswid=7467

In the article, the author states the problem with subs and eval. The way to get around it would be to use Maple as the symbolic engine but my goal is to not rely on Maple at all. Do you know of any sources that tackle this problem? How I would compute the function f and J, evaluate at a point, all numerically? Thank you for the link, I will read through it!

[Stephen23]

I think a core point is the randomness/arbitrariness of the functions. Presumably there are limits to their complexity, or perhaps you have some predetermined finite set of these functions? In both of these cases it would be possible to define the set functions once (in code) before starting this algorithm, and then simply perform a "random" selection from this set. This would allow symbolic functions, and yet would be much faster than generating the functions on the fly.

But is all depends on how random they need to be, and what knowledge you have in advance about these functions.

[Kristoffer Lindvall]

Yes exactly. The problem is I have to generate these functions and solve the system of equations before continuing to the next time interval. I don't know anything in advance.

[Stephen23]

"I don't know anything in advance." Not even the maximum order of the polynomials?

[Kristoffer Lindvall]

Of the solution polynomial yes (I construct it with weight coefficients), but I have no prior knowledge of the set of equations f and the jacobian J that I need to evaluate. From the solution of the equations f I get the weights for the solution polynomial.

[Kristoffer Lindvall]

As a side note, it is pretty interesting why Matlab's command subs is so slow. I don't understand why, considering that Maple runs it like no problem. I don't know all the differences between them, but Matlab should atleast be able to compete, but it doesn't.