How to encode a vector expression that depends on another vectorized expression?

Question

If, for example, I have three expressions: A , B and C as follows:

 A(i+1) = A(i) + C(i).k B(i+1) = B(i) + A(i).h C(i+1) = A(i) + B(i) 

where k and h are some constants, and m and n are the desired size of C i is the previous received value, i+1 is the next value. Now, if I use a for loop, then I can encode it as:

 A(1)= 2; B(1)= 5; C(1)= 3; for i=1:10 A(i+1) = A(i) + C(i)*2; B(i+1) = B(i) + A(i)*3; C(i+1) = A(i) + B(i); end 

And everything works fine. But I want to encode it in vector form , as well as without using a loop. But the problem is that I do not know how to get around the dependency:

• A by its previous value and previous value C
• B on it previous values ​​and previous values C A
• C for previous values ​​of A and B

Answer 1

Here, the matrix-based way to get the n value of the vector [A;B;C] . I would not call it vectorization, but this can significantly speed up the process:

 [A,B,C] = deal(zeros(11,1)); A(1)= 2; B(1)= 5; C(1)= 3; %% // Original method for k=1:10 A(k+1) = A(k) + C(k)*2; B(k+1) = B(k) + A(k)*3; C(k+1) = A(k) + B(k); end %% // Matrix method: %// [ A ] [1 0 2][ A ] %// | B | = |3 1 0|| B | %// [ C ] [1 1 0][ C ] %// i+1 i %// %// [ A ] [1 0 2][ A ] [1 0 2] ( [1 0 2][ A ] ) %// | B | = |3 1 0|| B | = |3 1 0| * ( |3 1 0|| B | ) %// [ C ] [1 1 0][ C ] [1 1 0] ( [1 1 0][ C ] ) %// i+2 i+1 i %// Thus, this coefficient matrix taken to the n-th power, multiplied by the input %// vector will yield the values of A(n+1), B(n+1), and C(n+1): M = [1 0 2 3 1 0 1 1 0]; isequal(M^10*[A(1);B(1);C(1)],[A(11);B(11);C(11)]) 

In fact, you can use M for the corresponding power (positive or negative) to get any [A,B,C] _n from any [A,B,C] _k ...

Answer 2

First, forgive me for abusing the Matlab syntax to express mathematical material.

Consider the following code, where we do the same as in your example. Note that A,B,C are strings of X

 X = zeros(3,N+1); X(:,1) = [2,5,3]; M= [1,0,2;3,1,0;1,1,0]; for i=1:N X(:,i+1) = M*X(:,i); end 

This is just a matrix vector designation of the above code. I think this is even slower. Note that we could also calculate: X(:,i+1) = M^i * X(:,1) , which is even slower.

Please note that we can use eigenvalue decomposition:

 [V,D] = eigs(M); X(:,i+1) = [V*D*inv(V)]^i * X; 

therefore

 X(:,i+1) = V*D^i*inv(V) * X; 

So V*D^i*inv(V) is an explicit formula for the i+1 1th term of X I suggest calculating them analytically and reconnecting the formula that you enter into your code.

EDIT: I wrote code that should be close to the analytical solution of the system, you can compare the runtime. It seems that ultimately preallocation with your first method is still the fastest if you need ALL TERMS. If you need only one of them, my proposed method will certainly be faster.

 clear;clc N = 10000000; tic A(1)= 2; B(1)= 5; C(1)= 3; A = zeros(1,N+1); B=A;C=A; for i=1:N A(i+1) = A(i) + C(i)*2; B(i+1) = B(i) + A(i)*3; C(i+1) = A(i) + B(i); end toc tic X = zeros(3,N+1); X(:,1) = [2,5,3]; M= [1,0,2;3,1,0;1,1,0]; for i=1:N X(:,i+1) = M*X(:,i); end toc tic M= [1,0,2;3,1,0;1,1,0]; [V,D]=eig(M); v=0:N; d=diag(D); B=bsxfun(@power,repmat(d,1,N+1),v); Y=bsxfun(@times,V * B, V \[2;5;3]); toc tic M= [1,0,2;3,1,0;1,1,0]; [V,D]=eig(M); v=0:N; d=diag(D); Y = ones(3,N+1); for i=1:N Y(:,i+1) = d.*Y(:,i); end Y=bsxfun(@times,V * B, V \[2;5;3]); toc