Discussion :

[Fleur du Pays]

Tout d'abord bonjour à tous, j'aimerais avoir votre aide concernant la programmation de la méthode d'Euler Implicite, je vous montre ce que j'ai effectué.

EDIT: j'ai supprimé mes programmes suite à la copie par certains élèves de ma promo (ils se reconnaitrons :p), demandez donc moi plutôt que de copier des programmes à moitié faux que vous ne comprenez pas ^^

[phryte]

Bonjour.

function y=dF(x);
y=[0 1;-1 0];

Tu n'utilises pas l'argument d'entrée x ?

[Fleur du Pays]

dF correspond à la Jacobienne de F, (sauf erreur de ma part sa donne celà ^^) donc je n'ai pas utilisé x pour cela, néanmoins je m'étais dit qu'en y mettant sous cette forme se serait plus simple si je devais le refaire avec une autre fonction.
Le fait que je ne l'utilise pas pourrait donc créer un soucis? =/

EDIT:

J'ai modifié une ligne de ma fonction de Newton par celle ci

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Nw=Nw-inv(feval(dF,Nw))*feval(F,Nw);

Cela m'enlève déjà la seconde erreur, mais la première persiste =/...

[Jerome Briot]

Je ne vois pas très bien l'intérêt de FEVAL ici

A mon avis, tu devrais reprendre ton algorithme depuis le début

[Fleur du Pays]

Je me suis inspiré justement de ce que j'ai pu trouver afin d'essayer de construire cette petite chose bancale qui ne marche pas ^^'... Comme je l'ai souligné, je suis plutôt débutant en la matière. Il me semble avoir réussi pour Euler Explicite, mais même en reprenant tout du début pour Euler Implicite, le problème est que je ne vois pas par où commencer différemment...

Concernant feval, je dois bien avouer ne pas trop connaitre les moments précis où sont utilisation est strictement nécessaire, mais les fonctions f et df sont des arguments d'entrée de la fonction newton et "il est donc
nécessaire d'utiliser la commande feval pour les appeler" à ce que j'ai pu recevoir comme informations jusqu'à présent...

Mis à part cela, il me semble par contre m'être rendu compte que je n'applique pas newton aux bonnes fonction F et dF pour construire un schéma d'Euler Implicite. Il me faudra relire quelque peu mon cours là dessus pour en être certain.

Mais j'ai quelque peu l'impression d'être dans une impasse,ne sachant plus trop que faire pour m'en sortir...

[salseropom]

Salut, dans ta méthode de Newton, tu peux pas plutôt faire plus simplement :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
Nw=Nw-F(Nw)\dF(Nw);

Ca m'a l'air plus simple.

Autre probème : il faut tester que ta matrice jacobienne soit inversible (donc 0 ne doit pas être une valeur propre...)
Si c'est le cas, tu peux t'inspirer d'une méthode de Newton modifiée... Si J(t) non inversible alors tu poses J(t) = J(t - 1) ou J(t) est la jacobienne à l'instant t

[FR119492]

Bonjour à tous.
salseropom a parfaitement raison: c'est non seulement plus simple, mais aussi plus rapide.
Pourquoi les gens ont-ils la manie d'inverser une matrice lorsqu'ils doivent résoudre un système linéaire?
Jean-Marc Blanc

[Fleur du Pays]

Pourquoi les gens ont-ils la manie d'inverser une matrice lorsqu'ils doivent résoudre un système linéaire?
Jean-Marc Blanc

Pourrais je savoir de quoi, ou plutôt à quelle ligne tu fais allusion?
Car depuis hier, selon les aides que je demande, bcp me répondent par ce genre de phrases, sans justement m'expliquer la chose, ce qui me permettrait justement d'avancer et de ne pas rester bloquer un WE entier là dessus. Lorsque les EC sont grévistes, il est malheureusement parfois difficile de s'autoformer...

[FR119492]

Salut!

Pourrais je savoir de quoi, ou plutôt à quelle ligne tu fais allusion?

Je comparais simplement la ligne de code que tu as donnée dans ton message d'hier à 18h25 et celle donnée par salseropom dans son message d'hier à 21h41. Pour les explications, voir mon tutoriel http://jmblanc.developpez.com/algori...mes-lineaires/
Jean-Marc Blanc

[Fleur du Pays]

J'ai finalement tout repris en m'aidant d'un autre code que j'avais trouvé sur internet et que je ne comprenais pas au départ, je l'ai modifié de manière à ce qu'il me paraisse (à mes yeux) plus simple.

Je poste ce que cela donne si jamais cela peut aider quelqu'un à l'avenir, il me semble avoir cette fois réussi à obtenir quelque chose qui fonctionne correctement.
Si vous avez néanmoins quelques remarques je suis preneur.
Dans le code d'origine il y avait par exemple "eye(lenght(x))" ) la place de la matrice identité que j'ai tappé directement dans le code, j'ai compris son utilité si ma fonction F était différente, mais ne connaissant au départ pa sla fonction "lenght", j'ai préféré opté pour cette solution.
J'avais également fait un erreur d'indices dans ma fonction "Energie" au départ.

[...]

[...]

En espérant que ce code ne soit pas vulgairement recopié par les autres éleves de ma promo pour mon projet ^^

[salseropom]

Pourrais je savoir de quoi, ou plutôt à quelle ligne tu fais allusion?
Car depuis hier, selon les aides que je demande, bcp me répondent par ce genre de phrases, sans justement m'expliquer la chose, ce qui me permettrait justement d'avancer et de ne pas rester bloquer un WE entier là dessus. Lorsque les EC sont grévistes, il est malheureusement parfois difficile de s'autoformer...

en gros, si tu veux calculer

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
x = A^(-1) * b

où A est une matrice carrée inversible, b un vecteur colonne et x l'inconnue, il vaut mieux chercher x tel que

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
A * x = b

d'où ce que disait Jean-Marc sur la résolution des systèmes linéaires.

D'un point de vue temps de calcul, calculer l'inverse d'une matrice prend un temps extrêmement long ... et en plus tu as plein d'erreurs numériques qui s'ajoutent (car tu as bcp d'opérations à faire).

A retenir : ne JAMAIS inverser une matrice mais TOUJOURS résoudre un système linéaire.

Selon moi, il y a un cas particulier : une matrice 2*2 (voire 3*3) où on peut calculer explicitement l'inverse d'une matrice. Mais cette situation est tellement rare...

[salseropom]

J'ai finalement tout repris en m'aidant d'un autre code que j'avais trouvé sur internet et que je ne comprenais pas au départ, je l'ai modifié de manière à ce qu'il me paraisse (à mes yeux) plus simple.

Je poste ce que cela donne si jamais cela peut aider quelqu'un à l'avenir, il me semble avoir cette fois réussi à obtenir quelque chose qui fonctionne correctement.
Si vous avez néanmoins quelques remarques je suis preneur.
Dans le code d'origine il y avait par exemple "eye(lenght(x))" ) la place de la matrice identité que j'ai tappé directement dans le code, j'ai compris son utilité si ma fonction F était différente, mais ne connaissant au départ pa sla fonction "lenght", j'ai préféré opté pour cette solution.
J'avais également fait un erreur d'indices dans ma fonction "Energie" au départ.





En espérant que ce code ne soit pas vulgairement recopié par les autres éleves de ma promo pour mon projet ^^

change ton

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
while(err > eps)

en

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
 
itermax = 100000;
iter = 1;
while( (err > eps) & (iter < maxiter) )
 
blabla;
iter = iter + 1;
end

il faut que tu mettes un cas de sortie en cas de divergence de ton algo

[FR119492]

Salut!

Selon moi, il y a un cas particulier

En fait, il y en a au moins un autre: si la matrice est orthogonale; dans ce cas, son inverse est identique à sa transposée.
Jean-Marc Blanc

[Fleur du Pays]

Merci pour l'amélioration de Newton. Je pense en avoir enfin terminé ^^

Une petite question tout de même, même si celle ci ne correspond plus au sujet initial, de quelle manière dois-je m'y prendre si je souhaite comparer les Méthodes de Runge Kutta et de Crank Nicholson (que j'ai déjà implémentés) pour la résolution d'une équa diff de la forme y"+sin(y)=0 (Pendule linéaire)?

J'avais pensé à tracer sur un même graphique ||Solution_schéma_CN(t) - Solution_exacte(t)|| et ||Solution_schéma_RK4(t) - Solution_exacte(t)||
Cela vous semble une méthode adéquate?

[salseropom]

Merci pour l'amélioration de Newton. Je pense en avoir enfin terminé ^^

Une petite question tout de même, même si celle ci ne correspond plus au sujet initial, de quelle manière dois-je m'y prendre si je souhaite comparer les Méthodes de Runge Kutta et de Crank Nicholson (que j'ai déjà implémentés) pour la résolution d'une équa diff de la forme y"+sin(y)=0 (Pendule linéaire)?

J'avais pensé à tracer sur un même graphique ||Solution_schéma_CN(t) - Solution_exacte(t)|| et ||Solution_schéma_RK4(t) - Solution_exacte(t)||
Cela vous semble une méthode adéquate?

Pour comparer des erreurs, ne prend pas l'erreur absolue (ce que tu fais) mais l'erreur relative (tant que la valeur exacte est différente de 0)

Autre point, tu peux aussi comparer les temps de calculs. Si tu prends une méthode d'ordre 4 à pas variable, tu feras plus de calculs qu'une bête méthode d'ordre 1

De mémoire, je crois que Crand Nicholson est une méthode d'ordre 1 (plus très sûr de moi) et runge kutta d'ordre ... tout dépend la méthode que tu choisis.

Sinon, tu peux aussi regarder du côté de ode45 et toutes ses amies

[Fleur du Pays]

Entendu, je vais essayer cela, pour simple info, j'ai opté pour du Runge Kutta d'ordre 4, et pour info il me semble que Crank nicholson est en fait d'ordre 2.
on m'avait sinon proposé une autre méthode que voici

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

[...]

Toute la première partie jusqu'au trait ne pose pas de pb, si j'ai bien compris nous allons comparer les résultats de chaque méthode avec la solution exacte en fonction de différents pas de temps (de 0.1 à 0.01 apparement)
Je l'ai d'ailleurs modifié par rapport à la version originale pour que cela fonctionne dans mon cas.
Mais c'est dans la dernière partie que je coince (plus précisément la partie en rouge), c'est la fonction polyfit qui me pose problème, car je ne la connais pas et n'ai pas réussi à bien saisir quel est son rôle. J'étasi malheureusement censé la voir en TP la semaine dernière, mais pour cause de greve cela n'a pas été...

Le résultat est alors le suivant: Qui permet bien de ce rendre compte que Runge Kutta d'ordre 4 est une meilleur approximation.

edit: Le graphique était faix, à cause d'une erreur dans ma solution exacte.
J'ai finalement compris l'interet du polyfit.

Merci à vous pour votre aide, je marque le tire du sujet comme résolu =)