aide de code MATLAB. Backward Euler méthode

Question

Voici le code MATLAB/FreeMat J'ai obtenu de résoudre un ODE numériquement en utilisant la méthode backward Euler. Cependant, les résultats sont incompatibles avec les résultats de mon manuel, et parfois même ridiculement incohérents. Quel est le problème avec le code?

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

Votre méthode est une méthode d'un nouveau type. Ce n'est ni l'arrière ni l'avant d'Euler.:-)

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

arrière Euler solve in y1: y1 - h*f(x1,y1) = y0

Votre méthode: y1 = y0 +h*f(x0,x0+h*f(x0,y0))

Votre méthode est pas Euler vers l'arrière.

Vous ne résolvez pas dans y1, vous venez d'estimer y1 avec la méthode Euler directe. Je ne veux pas poursuivre l'analyse de votre méthode, mais je crois qu'elle se comportera mal, même par rapport à Euler, puisque vous évaluez la fonction f au mauvais point.

Voici la méthode la plus proche de votre méthode que je puisse penser, explicite aussi, qui devrait donner de bien meilleurs résultats. Il est Heun's Method:

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

Answer 2

Le seul problème que je peux repérer est que la ligne:

n=(xfinal-xinit)/h 

devrait être:

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

Pour éviter un nombre de pas négatif. Si vous vous déplacez dans la direction négative-x, assurez-vous de donner à la fonction une taille de pas négative.

Vos réponses dévient probablement à cause de la façon dont vous grossièrement votre réponse d'approximation. Pour obtenir un résultat semi-précis, deltaX doit être très très petit et votre taille de pas doit être très très petite.

PS. Ce n'est pas la «méthode d'Euler rétrograde», c'est juste la méthode habituelle du vieux Euler.

Dans ce devoir s'il vous plaît marquer ainsi.

Answer 3

Jetez un oeil à numerical recipes, en particulier le chapitre 16, l'intégration des équations différentielles ordinaires. La méthode d'Euler est connu pour avoir des problèmes:

Il y a plusieurs raisons pour lesquelles la méthode d'Euler n'est pas recommandé pour une utilisation pratique, parmi eux, (i) la méthode n'est pas très précis par rapport à d'autres, plus fantaisistes, les méthodes fonctionnent à l'équivalent stepize, et (ii) n'est pas non plus très stable

Donc, à moins que vous sachiez que votre manuel utilise la méthode d'Euler, je ne m'attendrais pas à ce que les résultats correspondent. Même si c'est le cas, vous devez probablement utiliser une taille de pas identique pour obtenir un résultat identique.

Answer 4

À moins que vous vraiment voulez résoudre un ODE par la méthode d'Euler que vous avez écrit vous-même, vous devriez jeter un oeil à built-in ODE solvers.

Sur un sidenote: vous n'avez pas besoin de créer x(i) à l'intérieur de la boucle comme ceci: x(i+1) = x(i)+h;. Au lieu de cela, vous pouvez simplement écrire x = xinit:h:xfinal;. En outre, vous pouvez écrire n = round(xfinal-xinit)/h); pour éviter les avertissements.

---

Voici les solveurs implémentés par MATLAB.

ode45 est basé sur une formule explicite Runge-Kutta (4,5), la paire Dormand-Prince . Il est en une seule étape solveur - dans le calcul de y (tn), il faut seulement la solution à la suite précédentes point dans le temps, y (tn-1). En générale, ode45 est la meilleure fonction pour appliquer comme premier essai pour la plupart des problèmes .

L'ode23 est une implémentation d'une paire de Runge-Kutta explicite (2,3) de Bogacki et Shampine. Il peut être plus efficace que ode45 aux tolérances et en présence de rigidité modérée. Comme ode45, ode23 est un solveur en une étape.

L'ode113 est un solveur variable PECE Adams-Bashforth-Moulton. Il peut être plus efficace que ode45 à des tolérances strictes et lorsque la fonction de fichier ODE est particulièrement cher à évaluer. ode113 est un solveur multi-pas - il faut normalement les solutions à plusieurs points précédents pour calculer la solution actuelle .

Les algorithmes ci-dessus sont destinés à résoudre les systèmes non-rigides. Si elles apparaissent trop lentes, essayez l'une des solutions suivantes: .

ode15s est un solveur à ordre variable basé sur les formules de différentiation numérique (NDF). Optionnellement, il utilise les formules de différenciation vers l'arrière (BDF, également connu sous le nom de méthode Gear) qui sont généralement moins efficaces. Comme ode113, ode15s est un solveur multi-étapes. Essayez ode15s lorsque ode45 échoue, ou est très inefficace, et vous soupçonnez que le problème est raide, ou lors de la résolution un problème différentiel-algébrique.

ode23s est basé sur une formule Rosenbrock modifiée d'ordre 2. Comme il est un résolveur en une seule étape, il peut être plus efficace que ode15s à brut tolérances. Il peut résoudre certains types de problèmes rigides pour lesquels ode15s n'est pas efficace.

L'ode23t est une implémentation de la règle trapézoïdale utilisant un interpolant "libre" . Utilisez ce solveur si le problème n'est que modérément rigide et vous avez besoin d'une solution sans amortissement numérique . ode23t peut résoudre les DAE.

ode23tb est une implémentation de TR-BDF2, une implicite formule de Runge-Kutta avec une première étape qui est une étape de règle trapézoïdale et un second étage qui est une formule de différenciation arrière d'ordre deux. Par construction, la même matrice d'itération est utilisée pour évaluer les deux étages. Comme les odes23, ce solveur peut être plus efficace que les ode15 aux tolérances .

Answer 5

Je pense que ce code pourrait fonctionner. Essaye ça.

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

Le code est correct. Juste vous devez ajouter une autre boucle dans la boucle for. Pour vérifier le niveau de cohérence.

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 

J'ai vérifié une fonction fictive et les résultats étaient prometteurs.