un code Matlab aide. Procédé Backward Euler

StackOverflow https://stackoverflow.com/questions/2937183

Question

Voici le MATLAB® / Code de FreeMat je suis arrivé à résoudre un ODE numériquement en utilisant la méthode arrière Euler. Cependant, les résultats sont en contradiction avec les résultats de mes manuels scolaires, et parfois même ridiculement incompatibles. 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
Était-ce utile?

La solution

Votre méthode est une méthode d'un nouveau genre . Il est ni en arrière, ni Euler. : -)

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

Backward 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 en arrière.

Vous n'estimez pas résolu dans y1, vous venez y1 avec la méthode d'Euler avant. Je ne veux pas poursuivre l'analyse de votre méthode, mais je crois qu'il se comporte mal en effet, même par rapport à Euler, puisque vous évaluez la fonction f au mauvais.

Voici la méthode la plus proche de votre méthode que je peux penser, explicite aussi bien, ce qui devrait donner des résultats bien meilleurs. Il est Méthode de Heun :

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

Autres conseils

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 de négatif. Si vous vous déplacez dans la direction x négative, assurez-vous de donner la fonction d'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 la taille de votre étape doit être très très très faible.

PS. Ce n'est pas la « méthode d'Euler en arrière, » il est juste régulière ancienne méthode d'Euler.

Dans ce tag devoirs s'il vous plaît ainsi.

Jetez un oeil à recettes numériques , en particulier le chapitre 16, l'intégration des différentielles ordinaires équations. 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 ne sont pas très précis par rapport à d'autres, plus fantaisistes, les méthodes fonctionnent à l'StepSize équivalent, et (ii) ni est-il très stable

Donc, sauf si vous savez que votre manuel utilise la méthode d'Euler, je n'attendre les résultats de match. Même si elle est, vous avez probablement utiliser une taille de pas identique pour obtenir un résultat identique.

Sauf si vous vraiment voulez résoudre une EDO via la méthode d'Euler que vous avez écrit par vous-même, vous devriez jeter un oeil à intégré solveurs ODE .

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

Voici les solveurs mises en œuvre par Matlab.

  

ode45 est basé sur une explicite   formule de Runge-Kutta (4,5), la   paire Dormand-Prince. Il est une étape   solveur - dans le calcul de y (tn), il a besoin   seule la solution au immédiatement   point temporel précédent, y (tn-1). Dans   général, ode45 est la meilleure fonction   appliquer un premier essai pour la plupart   problèmes.

     

ode23 est une mise en œuvre d'un   explicite paire de Runge-Kutta (2,3) de   Bogacki et Shampine. Il peut être plus   efficace que ode45 au brut   tolérances et en présence de   raideur modérée. Comme ode45, ode23   est un résolveur en une seule étape.

     

ode113 est un ordre de grandeur   solveur Adams-Bashforth-Moulton PECE.   Il peut être plus efficace que ode45 à   tolérances strictes et lorsque le ODE   fonction de fichier est particulièrement   coûteux à évaluer. ode113 est un   solveur multi-étapes - normalement les besoins   les solutions à plusieurs précédentes   points de temps pour calculer le courant   Solution.

     

Les algorithmes ci-dessus sont destinés à   résoudre des systèmes non raides. Si elles apparaissent   être trop lent, essayez d'utiliser l'un des   les solveurs raides ci-dessous.

     

ode15s est un solveur de commande variable   basé sur la différenciation numérique   formules (NDF). En option, il utilise   les formules de différenciation vers l'arrière   (BDF, également connu sous le nom de la méthode de Gear)   qui sont généralement moins efficaces. Comme   ode113, ode15s est un solveur multi-étapes.   Essayez ode15s en cas d'échec ode45, ou est   très inefficace, et vous pensez que   le problème est raide, ou lors de la résolution   un problème d'écart algébrique.

     

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

     

ode23t est une mise en œuvre du   règle trapézoïdale en utilisant un « libre »   interpolée. Utilisez ce système de résolution si le   problème est que modérément raide et   vous avez besoin d'une solution sans numérique   amortissement. ode23t peut résoudre DAE.

     

ode23tb est une implémentation de   TR-BDF2, un Runge-Kutta implicite   formule avec un premier étage qui est un   l'étape de la règle trapézoïdale et une seconde   étape qui est un arrière   formule de différentiation d'ordre deux.   Par construction, la même itération   la matrice est utilisée en évaluant à la fois   étapes. Comme ode23s, ce solveur peut   être plus efficace que ode15s au brut   tolérances.

Je pense que ce code pourrait fonctionner. Essayez ceci. 

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

Le code est très bien. Juste vous devez ajouter une autre boucle dans la boucle. 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é pour une fonction factice et les résultats ont été prometteurs.

Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
scroll top