Come codificare un'espressione vettorizzata che dipende da altre espressioni vettorizzate?

Question

Se per esempio ho tre espressioni: A, B e C come segue:

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

dove k e h sono alcune costanti e m e n è la dimensione desiderata di C. i è il valore ottenuto in precedenza, i+1 è il valore successivo. Ora, se io uso for ciclo, quindi posso codice come:

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 

e funziona bene. Ma voglio codificarlo in un formato vettoriale , come in, senza dover utilizzare un ciclo. Ma il problema è che non so come aggirare la dipendenza di:

• A sul suo valore precedente e precedenti C valore
• B su di esso i valori precedenti e il valore precedente C di A
• C su i valori precedenti di A e B

Answer 1

in primo luogo, mi perdoni per aver abusato della sintassi Matlab per esprimere mathem roba atical.

Considera il seguente codice, dove facciamo esattamente lo stesso del tuo esempio. Si noti che A,B,C sono le righe di 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 

Questo è solo una notazione matrice vettore del codice precedente. Penso che sia ancora più lento. Nota che potremmo anche calcolare: X(:,i+1) = M^i * X(:,1) che è ancora più lento.

Si noti che possiamo usare la decomposizione autovalore:

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

Pertanto

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

Quindi V*D^i*inv(V) è una formula esplicita per il i+1 esimo termine di X. Suggerisco di calcolare quelle analiticamente e collegare di nuovo la formula che si ottiene nel codice.

EDIT: Ho scritto un codice che dovrebbe essere vicino alla soluzione analitica del sistema, è possibile confrontare i tempi di esecuzione. Sembra alla fine preallocazione con il tuo primo metodo è ancora il più veloce IF è necessario ALL i termini. Se hai solo bisogno di uno di loro, il mio metodo suggerito è sicuramente più veloce.

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 

Answer 2

Ecco un modo matriciale avere la n valore -esimo del vettore [A;B;C].Io non esattamente chiamare vettorializzazione, ma questo potrebbe accelerare le cose considerevolmente per voi:

[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 realtà è possibile utilizzare M alla potenza appropriata (positiva o negativa) per ottenere qualsiasi [A,B,C] _n da qualsiasi [A,B,C] _k ...