código MATLAB ayuda. método de Euler hacia atrás

Question

Aquí está el MATLAB / código que FreeMat llegaron a resolver un ODE numéricamente utilizando el href="http://en.wikipedia.org/wiki/Backward_Euler_method" rel="nofollow noreferrer"> hacia atrás . Sin embargo, los resultados son inconsistentes con los resultados de libros de texto, e incluso a veces ridículamente inconsistente. Lo que está mal con el código?

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

Su método es un método de un nuevo tipo . No es ni retrocede ni se transmita Euler. : -)

Reenviar Euler: y1 = y0 + h*f(x0,y0)

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

Su método: y1 = y0 +h*f(x0,x0+h*f(x0,y0))

Su método es no hacia atrás Euler.

no resuelven en y1, que acaba de estimar y1 con el método de Euler hacia adelante. No quiero seguir el análisis de su método, pero creo que va a comportarse mal hecho, incluso en comparación con Euler hacia adelante, ya que se evalúa la f función en el punto equivocado.

Este es el método más cercana a su método que se me ocurre, explícita, así, que debe dar resultados mucho mejores. Es Método de Heun :

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

Answer 2

El único problema que puedo observar es que la línea:

n=(xfinal-xinit)/h

En caso de ser:

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

Para evitar un recuento paso negativo. Si usted se está moviendo en la dirección negativa-x, asegúrese de dar la función de un tamaño negativo paso.

Sus respuestas probablemente se desvían por la forma en que se aproxima toscamente está su respuesta. Para obtener un resultado semi-precisa, deltaX tiene que ser muy muy pequeño y el tamaño de paso tiene que ser muy muy muy pequeña.

PS. Este no es el "método de Euler hacia atrás", es simplemente normal antiguo método de Euler.

Si esta es la tarea etiqueta Por favor, que así sea.

Answer 3

Para consultar numérica recetas , específicamente el capítulo 16, la integración de diferencial ordinaria ecuaciones. el método de Euler es conocida por tener problemas:

Hay varias razones por las que no se recomienda el método de Euler para el uso práctico, entre ellos, (i) el método no es muy preciso cuando se compara con otros, más elegantes, los métodos se ejecutan en el tamaño de paso equivalente, y (ii) tampoco es muy estable

Así que a menos que sepa que su libro de texto está usando el método de Euler, yo no esperaría los resultados de igualar. Incluso si lo es, es probable que tenga que utilizar un tamaño de paso idénticos para obtener un resultado idéntico.

Answer 4

A menos que realmente desea resolver una EDO a través del método de Euler que has escrito por sí mismo que debe tener una mirada en incorporado en solucionadores ODE .

En una nota: no es necesario para crear x(i) dentro del bucle de esta manera: x(i+1) = x(i)+h;. En lugar de ello, simplemente hay que escribir x = xinit:h:xfinal;. Además, es posible que desee escribir n = round(xfinal-xinit)/h); a las advertencias de evitar.

---

Aquí están los solucionadores implementadas por MATLAB.

ode45 se basa en una explícita Runge-Kutta (4,5) fórmula, el Dormand-Prince par. Es un solo paso solucionador - en el cálculo y (tn), que necesita sólo la solución en el inmediatamente punto precedente tiempo, y (tn-1). En en general, ode45 es la mejor función para aplicar como un primer intento para la mayoría problemas.

ode23 es una implementación de una explícita de Runge-Kutta (2,3) par de Bogacki y Shampine. Puede ser más eficiente que ode45 en crudo tolerancias y en presencia de rigidez moderada. Al igual que ode45, ode23 es un programa de solución de un solo paso.

ode113 es un orden variable de solucionador de Adams-Moulton Bashforth-PECE. Puede ser más eficiente que en ode45 tolerancias estrictas y cuando la ODE función de archivo es particularmente costosa de evaluar. ode113 es una MultiStep solucionador - Normalmente necesidades las soluciones a varios precedentes puntos de tiempo para calcular la corriente solución.

Los algoritmos anteriores tienen la intención de resolver sistemas nonstiff. Si aparecen a ser excesivamente lenta, intente utilizar una de los solucionadores rígidos a continuación.

ode15s es un solucionador de orden variable de basado en la diferenciación numérica fórmulas (NDFs). Opcionalmente, se utiliza las fórmulas de diferenciación atrasadas (BDFs, también conocido como el método de Gear) que son por lo general menos eficiente. Me gusta ode113, ode15s es un solucionador de varios pasos. Trate ode15s cuando ode45 falla, o es muy ineficiente, y se sospecha que el problema es dura, o en la resolución de un problema diferencial-algebraica.

ode23s se basa en una modificación fórmula Rosenbrock de orden 2. Debido es un programa de solución de una sola etapa, puede ser más eficiente que ode15s en crudo tolerancias. Puede resolver algunos tipos de problemas de rigidez para los que no es ode15s eficaz.

ode23t es una implementación de la regla trapezoidal utilizando un "libre" interpolador. Utilice este solucionador de si el problema es sólo moderadamente rígido y que necesita una solución sin numérica mojadura. ode23t puede resolver DAE.

ode23tb es una implementación de TR-BDF2, una implícita de Runge-Kutta fórmula con una primera etapa que es una paso regla trapezoidal y una segunda etapa que es un retroceso fórmula diferenciación de orden dos. Por construcción, la misma iteración matriz se utiliza en la evaluación de tanto stages. Al igual que ode23s, este solucionador de mayo ser más eficiente que ode15s en crudo tolerancias.

Answer 5

creo que este código podría funcionar. Prueba esto.

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

El código está muy bien. Sólo hay que añadir otro bucle dentro del bucle. Para comprobar el nivel de consistencia.

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

I comprueba para una función maniquí y los resultados fueron prometedores.