MATLAB help codice. metodo di Eulero all'indietro

Question

Questa è la MATLAB / FreeMat codice che ho avuto modo di risolvere un ODE numericamente usando l'href="http://en.wikipedia.org/wiki/Backward_Euler_method" rel="nofollow noreferrer"> all'indietro metodo . Tuttavia, i risultati sono in contrasto con i miei risultati di libri di testo, e talvolta anche ridicolmente inconsistenti. Cosa c'è di sbagliato con il codice?

function [x,y] = backEuler(f,xinit,yinit,xfinal,h)

    %f - this is your y prime
    %xinit - initial X
    %yinit - initial Y
    %xfinal - final X
    %h - step size

    n = (xfinal-xinit)/h; %Calculate steps

    %Inititialize arrays...
    %The first elements take xinit and yinit corespondigly, the rest fill with 0s.
    x = [xinit zeros(1,n)];
    y = [yinit zeros(1,n)];

    %Numeric routine
    for i = 1:n
        x(i+1) = x(i)+h;
        ynew = y(i)+h*(f(x(i),y(i)));
        y(i+1) = y(i)+h*f(x(i+1),ynew);
    end
end

Answer 1

Il metodo è un metodo di di un nuovo tipo . Non è né all'indietro né trasmettere Eulero. : -)

Inoltra Eulero: y1 = y0 + h*f(x0,y0)

Eulero solve in y1: y1 - h*f(x1,y1) = y0

Il tuo metodo: y1 = y0 +h*f(x0,x0+h*f(x0,y0))

Il tuo metodo è non di Eulero all'indietro.

Non c'è risolvere in y1, basta stimare y1 con il metodo di Eulero in avanti. Non voglio proseguire l'analisi del tuo metodo, ma credo che si comporterà male davvero, anche a fronte di Eulero in avanti, dal momento che il calcolo della funzione f nel punto sbagliato.

Ecco il metodo più vicino alla tua metodo che mi viene in mente, esplicita così, che dovrebbe dare risultati molto migliori. E ' Metodo di Heun :

y1 = y0 + h/2*(f(x0,y0) + f(x1,x0+h*f(x0,y0)))

Answer 2

L'unico problema che può individuare è che la linea:

n=(xfinal-xinit)/h

Dovrebbe essere:

n = abs((xfinal-xinit)/h)

Per evitare un conteggio dei passi negativo. Se ci si muove in direzione negativa-x, assicurarsi di dare alla funzione una dimensione negativa passo.

Le risposte probabilmente si discostano a causa di come grossolanamente si sta approssimando la vostra risposta. Per ottenere un risultato semi-accurata, deltaX deve essere molto molto piccolo e le dimensioni passo deve essere molto molto molto piccolo.

PS. Questo non è il "metodo di Eulero," si tratta solo di regolare vecchio metodo di Eulero.

Se questo è compito tag prega così.

Answer 3

Date un'occhiata al numerica ricette , in particolare il capitolo 16, l'integrazione dei differenziali ordinarie equazioni. il metodo di Eulero è noto per avere problemi:

  

Ci sono diversi motivi che il metodo di Eulero non è raccomandato per l'uso pratico, fra loro, (i) il metodo non è molto preciso rispetto ad altri, più elaborate, metodi funzionare alla stepSize equivalente, e (ii) non è nemmeno molto stabile

Quindi, se non si conosce il vostro libro di testo sta usando il metodo di Eulero, non vorrei aspettare i risultati a partita. Anche se si tratta, probabilmente dovete usare una dimensione passo identico per ottenere un risultato identico.

Answer 4

A meno che davvero vuole risolvere un ODE tramite il metodo di Eulero che hai scritto da solo si dovrebbe avere uno sguardo a built-in risolutori ODE .

Su un sidenote: non è necessario creare x(i) all'interno del ciclo in questo modo: x(i+1) = x(i)+h;. Invece, si può semplicemente scrivere x = xinit:h:xfinal;. Inoltre, si consiglia di scrivere n = round(xfinal-xinit)/h); agli avvertimenti evitare.

---

Ecco le risolutori attuati da MATLAB.

  

ode45 è basato su un esplicito   Runge-Kutta (4,5) formula, la   pair Dormand-Prince. Si tratta di un one-step   risolutore - nella computazione y (tn), è necessario   solo la soluzione al immediatamente   precedente punto di tempo, y (tn-1). Nel   generale, ode45 è la migliore funzione di   applicare come primo tentativo per la maggior parte   i problemi.

     

ode23 è un'implementazione di un   esplicito Runge-Kutta (2,3) coppia di   Bogacki e Shampine. Si può essere più   efficiente di ode45 al greggio   tolleranze e in presenza di   rigidità moderata. Come ode45, ode23   è un risolutore a singolo passo.

     

ode113 è un ordine variabile   Adams-Bashforth-Moulton PECE risolutore.   Può essere più efficiente di ode45 a   tolleranze severe e quando l'ODE   funzione di file è particolarmente   costoso da valutare. ode113 è un   Multistep risolutore - normalmente esigenze   le soluzioni a vari precedenti   punti di tempo per calcolare la corrente   soluzione.

     

Gli algoritmi di cui sopra sono destinate a   risolvere sistemi nonstiff. Se appaiono   ad essere eccessivamente lenta, provare a utilizzare uno dei   i risolutori rigidi sotto.

     

ode15s è un risolutore di ordine variabile   basato sulla differenziazione numerica   formule (NDF). Opzionalmente, utilizza   le formule di differenziazione all'indietro   (BDFs, noto anche come il metodo di Gear)   che di solito sono meno efficienti. Piace   ode113, ode15s è un risolutore multistep.   Prova ode15s quando ode45 non riesce, o è   molto inefficiente, e si sospetta che   il problema è rigida, o quando la soluzione   un problema differenziale-algebrico.

     

ode23s si basa su una versione modificata   Rosenbrock formula di ordine 2. Poiché   è un risolutore a singolo passo, può essere   più efficiente rispetto a ode15s greggio   tolleranze. E 'in grado di risolvere alcuni tipi di   problemi rigidi per i quali non è ode15s   efficace.

     

ode23t è un'implementazione del   regola trapezoidale utilizzando un "libero"   interpolante. Utilizzare questa risolutore se il   problema è solo moderatamente rigido e   avete bisogno di una soluzione senza numerica   smorzamento. ode23t può risolvere DAE.

     

ode23tb è un'implementazione del   TR-BDF2, un implicito Runge-Kutta   formula con un primo stadio che è un   regola di passaggio trapezoidale ed una seconda   fase che è arretrato   differenziazione formula di ordine due.   Per costruzione, la stessa iterazione   matrice viene utilizzato per valutare sia   stages. Come ode23s, questo risolutore può   essere più efficiente rispetto a ode15s greggio   tolleranze.

Answer 5

Credo che questo codice potrebbe funzionare. Prova questo.

for i =1:n
    t(i +1)=t(i )+dt;
    y(i+1)=solve('y(i+1)=y(i)+dt*f(t(i+1),y(i+1)');
    end 

Answer 6

Il codice va bene. Solo si deve aggiungere un altro ciclo all'interno del ciclo for. Per controllare il livello di coerenza.

if abs((y(i+1) - ynew)/ynew) > 0.0000000001 
    ynew = y(i+1);
    y(i+1) = y(i)+h*f(x(i+1),ynew);
end

ho controllato per una funzione fittizia ed i risultati sono stati promettenti.