2015-02-24 21 views
7

ho la seguente equazione che voglio risolvere rispetto a:Matlab: Soluzione un'equazione logaritmica

x = (a-b-c+d)/log((a-b)/(c-d)) 

dove x, sono noti b, c e d. Ho usato Wolfram Alpha per risolvere l'equazione, e il risultato è:

a = b-x*W(-((c-d)*exp(d/x-c/x))/x) 

dove W è la è la funzione log prodotto (Funzione W di Lambert). Potrebbe essere più facile vederlo allo Wolfram Alpha page.

Ho utilizzato la funzione incorporata di Matlab lambertW per risolvere l'equazione. Questo è piuttosto lento, ed è il collo di bottiglia nel mio script. C'è un altro, più veloce, modo di fare questo? Non deve essere preciso fino al decimo decimale.

MODIFICA: Non avevo idea che questa equazione sia così difficile da risolvere. Ecco un'immagine che illustra il mio problema. Le temperature b-d più LMTD variano in ogni fase temporale, ma sono note. Il calore viene trasferito dalla linea rossa (CO2) alla linea blu (acqua). Devo trovare la temperatura "a". Non sapevo che fosse così difficile da calcolare! : P enter image description here

+0

Con quale frequenza e in quale contesto è necessario risolvere questa equazione per 'a'? Come parte di qualche altro risolutore? I valori di 'x, b, c, ...' che vuoi valutare sono già noti? – knedlsepp

+0

Non so nulla della funzione lambert W. Tuttavia, non ho trovato la tua risposta molto utile, IrrationalPerson. @ knedlsepp: Sto usando per risolvere uno scambiatore di calore in una pompa di calore. Lo sto usando in qualche altro risolutore. Immagino che la funzione sia utilizzata almeno 30.000 volte, molto probabilmente più (sto facendo calcoli orari per un anno). Come menzionato, x, b, c e d sono noti. – ROLF

+0

Potresti modificare il tuo post per fornire alcuni valori realistici per 'x, b, c, d'? Il mio primo pensiero ingenuo è di provare a usare 'fzero()' per trovare la/e radice/i di '(ab-c + d)/ln ((ab)/(cd)) - x', ma questo può portare a risposte complesse per alcuni dei valori dei parametri casuali che ho provato. – eigenchris

risposta

4

Un'altra opzione è basata sulla semplice Wright ω function:

a = b - x.*wrightOmega(log(-(c-d)./x) - (c-d)./x); 

purché d ~= c + x.*wrightOmega(log(-(c-d)./x) - (c-d)./x) (cioè, d ~= c+b-a, x è 0/0 in questo caso). Questo è equivalente al ramo principale di Lambert W function, W, che penso sia il ramo di soluzione desiderato.

Proprio come con lambertW, c'è una funzione wrightOmega nella casella degli strumenti di matematica simbolica. Sfortunatamente, questo probabilmente sarà anche lento per un gran numero di input. Tuttavia, è possibile utilizzare il mio wrightOmegaq su GitHub per ingressi a virgola mobile a valori complessi (doppio o singolo). La funzione è più accurata, completamente vettorializzata e può essere da tre a quattro ordini di grandezza più veloce rispetto all'utilizzo dello wrightOmega integrato per gli ingressi a virgola mobile.

Per chi fosse interessato, wrightOmegaq si basa su questo eccellente carta:

Piers W. Lawrence, Robert M. Corless, e David J. Jeffrey, "Algorithm 917: Complex Double-Precision Evaluation of the Wright omega Function," Rapporti ACM su Software matematico, vol. 38, n. 3, articolo 20, pagg. 1-17, aprile 2012.

Questo algoritmo va oltre la convergenza cubica del Halley's method utilizzato in Cleve Moler di Lambert_W e utilizza un metodo radice conoscitiva quarto ordine convergenza (Fritsch, Shafer, & Crowley, 1973) a convergere in non più di due iterazioni.

Inoltre, per velocizzare ulteriormente il numero Lambert_W di Moler utilizzando espansioni in serie, vedere my answer at Math.StackExchange.

+0

Ho eseguito alcuni test sulle funzioni per ottenere una panoramica migliore: nell'intervallo '[-1/e, 0]' i tuoi valori iniziali potrebbero accelerare il processo, dato che l'implementazione di * Moler * richiede circa 6 iterazioni (al contrario a due). In pratica però è già abbastanza veloce, in quanto è completamente vettorializzato.Una cosa della tua risposta: non so se l'algoritmo di FSC che hai citato sia di ordine ancora più alto, ma il metodo di Halley è un metodo del terzo ordine. – knedlsepp

+1

@knedlsepp: Il metodo di Halley è un * secondo ordine * [metodo classe proprietario di casa] (http://en.wikipedia.org/wiki/Householder%27s_method). Il tasso di convergenza è superiore a quello dell'ordine, quindi il metodo di Halley ha * convergenza cubica *. La funzione 'Lambert_W' richiederà più iterazioni per raggiungere una determinata tolleranza per determinati valori critici. Il codice di esempio nella mia risposta Math.StackExchange è solo un'illustrazione dell'uso di espansioni in serie, non del codice da utilizzare per confrontare la velocità con un'altra funzione. La trama mostra che l'approccio più ingenuo usa più iterazioni, in particolare attorno a '-1/exp (1)'. – horchler

+0

Oh. * il metodo della classe Householder del secondo ordine: metodo di Halley *? Ciò potrebbe essere fonte di confusione, poiché solitamente un metodo * nth-order * riguarda l'ordine di convergenza. (Anche il documento FSC non sta utilizzando direttamente il metodo Householder applicato a questo problema, ma fornisce due algoritmi di respingente terzo e quarto ordine (di convergenza). – knedlsepp

1

Due (combinabili) opzioni:

  • È il vostro copione già Vectorized? Valutare la funzione per più di un singolo argomento. L'esecuzione di for i = 1:100, a(i)=lambertw(rhs(i)); end è più lenta di a=lambertw(rhs).
  • Se avete a che fare con il ramo a valori reali di LambertW (vale a dire i tuoi argomenti sono nell'intervallo [-1/e, inf)), è possibile utilizzare l'applicazione delle Lambert_W presentato dal Cleve Moler sul File Exchange.
+1

Buona presa sulla funzione 'Lambert_W' di Cleve Moler. Stavo per raccomandare una combinazione di utilizzo di 'exp' e' fzero', ma Moler sta facendo essenzialmente la stessa cosa ma con un metodo cubico, quindi è più che probabile più veloce. – TroyHaskin

+0

@TroyHaskin: Yup, ha dato una rapida occhiata al suo post sul blog in matematica. È abbastanza insolito per me vedere i metodi di ordine superiore rispetto all'iterazione di Newton. – knedlsepp

+0

@knedlsepp: Grazie per le risposte. Si prega di vedere la mia domanda aggiornata. – ROLF

Problemi correlati