[Matlab] Méthode de Newton-Raphson Résolu/Fermé

Question

Bonjour,

J'ai besoin de votre aide pour trouver les racines d'une équation:
 j'ai deux équations Ia = f(Ua) et Ib = f(Ub), je veux trouver les valeurs de Ua et Ub de tel façon que Ia-Ib = 0 ou d'une autre manière Ia = Ib, dans l'article est conseillé d'utiliser la méthode de newton raphson 
Mais en cherchant, tjrs sur internet, des infos sur cette méthode, j'ai pas compris grand chose.

Je peux vous donner l'expression de Ia et Ib en fonction de Ua et Ub respectivement.


Est ce que quelqu'un pourrait m'aider?.
Merci beaucoup d'avance pour votre aide

Sacabouffe · Accepted Answer

Salut attifa

Ben oui, si c'est pas ton domaine, c'est normal que ça te prenne un peu de temps pour comprendre ;-)
C'est comme quand j'ai lu ton code, bon, je voyais à peu près ce qu'il faisait mais ça ma parlait pas plus que ça ; l'électrophysiologie je connais pas et si veux essayer de voir à peu près ce qui se passe physiquement dans tout ce que t'as écrit, je risque d'y passer des heures :-D

Bonne nuit

Sacabouffe · Answer

Salut

Ben déjà Ua=Ub me paraît être une solution raisonnable :-D

Mais sinon, la méthode Newton-Raphson est décrite de manière plutôt claire sur Wikipédia :
https://fr.wikipedia.org/wiki/M%C3%A9thode_de_Newton#La_m.C3.A9thode

Mais t'as pas besoin d'implémenter ça, utilise les fonctions prédéfinies de Matlab, y a fzero par exemple. Tape help fzero ou doc fzero pour savoir comment elle s'utilise.

Ciao

atiffa · Answer

Bonjour sacabouffe:

Merci pour ta réponse, mais je pense que la fzero c'est pour résoudre un équation avec une seule variable (dans mon cas il y a Ua et Ub)
pour la méthode de newton raphson j'ai vraiment du mal a la comprendre (je ne suis pas une vraie matheuse
 :-(

C'est pas grave, on tous cas je te remercie infiniment pour ta disponibilité et pour ton aide

Cordialement
Attifa

Sacabouffe · Answer

Salut attifa

Ah oui pardon, j'avais pas pensé à ça :-(
Pour les fonctions de plusieurs variables, il y a fsolve, mais elle fait partie de l'Optimization Toolbox, je sais pas si tu l'as.

Cela dit, je sais pas si tu trouveras quelque chose pour résoudre simplement ton problème. En effet, en fait tu dois  trouver les solutions d'une équation de ce type g(x,y)=0.

Or si ta fonction g est une fonction continue admettant une dérivée partielle par rapport à y continue, d'après le théorème des fonctions implicites, si (x0,y0) est une solution de g(x0,y0)=0 telle que &#8706;g(x0,y0)/&#8706;y0&#8800;0, alors il existe un voisinage de U de x0, un voisinage V de y0 et une fonction continue &#966;:U&#8594;V tels que &#8704; x &#8712; U, g(x,&#966;(x))=0.

Pour simplifier dans ton cas particulier, comme la fonction g dont tu dois trouver les zéros est en fait définie par g(x,y)=f(x)-f(y), si ta fonction f est de classe C&#185;, pour toute solution (x0,y0) de f(x0)=f(y0) telle que f'(y0)&#8800;0, il existe des intervalles ]a,b[, ]c,d[ tels que x0 &#8712; ]a,b[, y0 &#8712; ]c,d[ et une fonction continue &#966;:]a,b[&#8594;]c,d[ tels que &#8704; x &#8712; ]a,b[, f(x)=f(&#966;(x)).

Tout ça pour dire qu'à moins d'un miracle (dérivée nulle à chaque fois, et encore...), les solutions de ton équation ne sont pas isolées ; pour chaque solution, il existe un voisinage de cette solution dans lequel il y a une infinité de solutions. 

Bonne journée

Char Snipeur · Answer

disons que la fonction à utiliser dans ton cas, c'est:
g(x,y)=f(x)-f(y)
comme toi, je ne penses pas que fzero soit adapter à cette résolution à deux variables. Il faut donc utiliser la méthode donné dans "généralisation" de la page Wikipedia.
Dans ton cas, le jacobien est assez simple. Cette méthode n'est pas très compliqué, il suffit de passer un peu de temps pour l'appliquer.

Sacabouffe · Answer

C'est pas un système carré d'équations non linéaires :-(
Une seule équation et deux inconnues.
Je pense pas que ça donnera quelque chose.

J'ai tenté l'expérience avec fslove (à la place de fzero) et une fonction test, Matlab utilise l'algorithme de Levenberg-Marquardt.

Sans l'Optimization Toolbox, y a possibilité d'utiliser fminsearch sinon et de minimiser |g|, mais l'éventualité de tomber sur un minimiseur local (x,y) de g tel que |g(x,y)|&#8800;0 (et donc f(x)&#8800;f(y)) est pas à exclure, bref... la misère de toute façon... :-D

Mais quelle que soit la méthode, en partant d'un initial guess donné, je sais pas trop quelle solution va sortir :-(
La solution minimisant la distance entre l'initial guess et la courbe représentative de la fonction implicite associée à la solution ?

atiffa · Answer

Voila sacabouffe,

Pour voir un peu a quoi ressemble mes équations, c'est sa mon code,

%%%%%%%%%%%%%%%%%%%%%%%

global cnaa cka ccla cx1a cx2a cnab ckb cclb cx1b cx2b 
global pnaa pka pnab pkb pclb pnaj pkj pclj   
global lpa lpb Aa Ab Aj RT faraday 
global zna zk zcl gammaMem


%Initial concentrations:

cnac0 = 7;     % mM
ckc0 = 72;     % mM
cclc0 = 16;    % mM
cx1c0 = 63;    % mM
cx2c0 = 142;   % mM

cnaa = 104;    % mM
cka = 5.3;     % mM
ccla = 102;    % mM
cx1a = 7.3;    % mM
cx2a = 81.4;   % mM

cnab = 104;    % mM
ckb = 5.3;     % mM
cclb = 102;    % mM
cx1b = 7.3;    % mM
cx2b = 81.4;   % mM

%Initial potentiel:

Ea0 = 2.5;
Eb0 = -3;
Et0 = Ea0 + Eb0;

%Permeabilties:

pnaa = 100e-9;       % cm./sec
pka = 50e-9;         % cm./sec

pnab = 20e-9;        % cm./sec
pkb = 463e-9;        % cm./sec
pclb = 541e-9;       % cm./sec

pnaj = 3e-9;         % cm./sec
pkj = 3e-9;          % cm./sec
pclj = 3e-9;         % cm./sec

%Hydraulic conductances:

lpa = 1e-12;  %cm3./dyne s  %cm3./(g*cm./s^2)*s %cm3./(1e-5N)*s
lpb = 1e-11;  %cm3./dyne s  %cm3./(g*cm./s^2)*s %cm3./(1e-5N)*s

%geometrical parameters:

Aa = 1.8;     %cm2/cm2 tissue
Ab = 8.8;     %cm2/cm2 tissue
Aj = 0.001;   %cm2/cm2 tissue, Alan Weinstein 2005
Vc0 = 0.001;  %cm3/cm2 tissue

%Parametrs

RT = 19.36;                      % R : the gas constant; T:the absolute temperature(RT: mmHg.l/mmol)
faraday = 96.48;                 % Faraday's constant (coul/mmol = mcoul/umol1.0)
gammaMem =1e-6;                  % memCapacitance = 1e-6 farads/cm2
% Les valences:

zna = 1;
zk = 1;
zcl = -1;


tspan = 0:1000;

y0 = [cnac0; ckc0; cclc0; cx1c0; cx2c0; Ea0; Eb0; Vc0];
[T,Y] = ode15s(@diffs, tspan, y0);

 
 cnac = Y(:,1);
 ckc = Y(:,2);
 cclc = Y(:,3);
 cx1c = Y(:,4);
 cx2c = Y(:,5);    
 Ua = Y(:,6)*faraday./RT;
 Ub = Y(:,7)*faraday./RT;
 Vc = Y(:,8);  
 
 subplot(3,1,1);plot(T,ckc,'b');
 subplot(3,1,2);plot(T,cclc,'r');
 subplot(3,1,3);plot(T,cnac,'g');

%%%%%%%%%%%%%%%%%

function dydt = diffs(t,y)

global cnaa cka ccla cx1a cx2a cnab ckb cclb cx1b cx2b 
global pnaa pka pnab pkb pclb pnaj pkj pclj   
global lpa lpb Aa Ab Aj RT faraday 
global zna zk zcl gammaMem

 % y = [cnac; ckc; cclc; cx1c; cx2c; Ea; Eb; Vc];
 
 cnac = y(1);
 ckc = y(2);
 cclc = y(3);
 cx1c = y(4);
 cx2c = y(5);    
 Ua = y(6)*faraday./RT;
 Ub = y(7)*faraday./RT;
 Vc = y(8);
 
 Ut = Ua + Ub;
  
 
 %%%l'expression de la concentration du chlore à partir des équations d'électroneutralités
 ccla = cnaa + cka - cx1a ;

 %%%Le flux ionique couplés
 [jnap jkp] = pumpnak(cnac,ckb);
 
 %%%les equations d'osmolalite:
 osmoa = cnaa + cka + ccla +  cx1a + cx2a;
 osmoc = cnac + ckc + cclc +  cx1c + cx2c;
 osmob = cnab + ckb + cclb +  cx1b + cx2b;
 
 %%%les equations de flux d'eau
 jwa = (lpa*RT).*(osmoc - osmoa);
 jwb = (lpb*RT).*(osmoc - osmob);
 
 
 %%% le flux des ions à travers la membrane apicale de la cellule du TCD
 jnaa = Aa*GHKflux(pnaa,zna,Ua,cnaa,cnac);
 jka = Aa*GHKflux(pka,zk,Ua,cka,ckc);
 
 %%% le flux des ions à travers la membrane basale de la cellule du TCD 
 jnab = Ab*GHKflux(pnab,zna,Ub,cnac,cnab) + jnap;
 jkb = Ab*GHKflux(pkb,zk,Ub,ckc,ckb) + jkp;
 jclb = Ab*GHKflux(pclb,zcl,Ub,cclc,cclb);
 
 %%% le flux des ions à travers l'espace transcellulaire
 jnaj = Aj*GHKflux(pnaj,zna,Ut,cnaa,cnab);
 jkj = Aj*GHKflux(pkj,zk,Ut,cka,ckb);
 jclj = Aj*GHKflux(pclj,zcl,Ut,ccla,cclb);
 
 
 %%% conservation de l'eau et des solutes dans le compartiment c
 dvcdt = Aa*jwa + Ab*jwb;
 dcnacdt =(Aa*jnaa - Ab*jnab - cnac.*(Aa*jwa + Ab*jwb))./Vc; 
 dckcdt =(Aa*jka - Ab*jkb - ckc.*(Aa*jwa + Ab*jwb))./Vc; 
 dcclcdt =(- Ab*jclb - cclc.*(Aa*jwa + Ab*jwb))./Vc; 
 dcx1cdt = (- cx1c.*(Aa*jwa + Ab*jwb))./Vc; 
 dcx2cdt = (- cx2c.*(Aa*jwa + Ab*jwb))./Vc; 

 Ia = faraday*Aa*(zna*(jnaa+jnaj) + zcl*jclj + zk*(jka+jkj));
 Ib = faraday*Ab*(zna*(jnab+jnaj) + zcl*(jclb+jclj) + zk*(jkb+jclj));
 
 dEadt = Ia/(gammaMem*Aa);
 dEbdt = Ib/(gammaMem*Ab);
 
 dydt = [ dvcdt;
          dcnacdt ; 
          dckcdt ; 
          dcclcdt ; 
          dcx1cdt ; 
          dcx2cdt ; 
          dEadt ;
          dEbdt ]; 
end
%%%%%%%%%%%%%%%%%%%%%%%%


%Input: cnac: The concentration of intracellular Na+ (mM)
%       ckb: The concentraton of serosal K+ (mM)

%Output: jnap: flux attributable to the sodium pump(Na+ efflux)(umol/cm2.sec)
%        jkp : flux attributable to the sodium pump(K+ influx) (umol/cm2.sec)

%References: Richard Latta and all.General Method for the Derivation and Numerical Solution 
%            of Epithelial Transport Models.J. Membrane Biol. 82, 67-82 (1984) 
 
F = 96.48;                 % Faraday's constant (coul/mmol = mcoul/umol1.0)
Imax = 5.61e-3;            % The maximum observable pump current(nA/cm2)
Kna = 14.2;                % The concentration of intracellular Na+ producing half-maximal pump current(mM)
Kk = 2.3;                  % The concentration of serosal K+ producing half-maximal pump current(mM)


Ip = Imax./((1+(Kna./cnac).^3).*(1+(Kk./ckb).^2));

jnap = 3*Ip/F;
jkp = - 2*Ip/F;

end
%%%%%%%%%%%%%%%%%%%%%%

function Jghkipq= GHKflux(pipq,zi,Upq,cip,ciq)

% jghk the Goldman - Hodgkin - Kartz equation: which was applied to all the
% conductive pathways (ion channels and paracellular pathways)
%
% Input: pipq: permeability in cm/s = cm3/(cm2.sec)
%        zi: valence of i solute
%        Upq: normalized membrane potential of compartment p and q (unitless)
%        cip,ciq: concentration of solute i in compartment p and q(mM)
%        
% Output: Jghkipq: flux of solute i from compartment p to compartment q (umol/cm2.s)
% Note: i = Na+,K+,Cl-
%       p,q = mucosal, cellular and cerosal compartment
% Reference: Richard Latta and all, 
%            J Membrane Biol 82,67-82,1984


 if (abs (Upq)<0.001)
     Jghkipq = pipq*(cip - ciq);
 else
    Jghkipq = pipq*zi*Upq.*((cip - ciq).*exp(-zi.*Upq))./(1-exp(-zi.*Upq));
 end
 
 end 
 
%%%%%%%%%%%%%%%%%%%%%%

J'ai commencer à avoir des sorties graphiques mais il faut ajuster les parametres pour obtenir les vrais courbes

J'ai besoin de varier les valeurs de Ua et Ub pour que a t=0 : Ia = Ib

Merci sacabouffe et merci Char Snipeur pour votre aide

Cordialement
Atiffa

Sacabouffe · Answer

Salut
Vois avec les infos du message 5, à la seule différence que le problème est pas tout à fait celui que tu présentais dans ton sujet.
Il te faut les solutions d'une équation du type f(x)=g(y).
Tu peux aussi essayer un logiciel de calcul formel pour voir essayer d'exprimer la relation explicite entre x et y, c'est-à-dire avec un truc du genre Maple, MuPAD, ou peut-être même la Symbolic Math Toolbox de Matlab, faut voir...
Ciao

atiffa · Answer

Finalement sacabouffe,

Je suis coincé, je ne sais plus par où commencer (cela m'angoisser vraiment...)

Merci pour ton aide c'est gentil de ta part

atiffa · Answer

Oui avant que j'oublie est ce que t'auras pas des code matlab avec la méthode de newton raphson (j'ai des mal à la comprend)


Merci d'avance

Sacabouffe · Answer

J'avais pas pensé, j'aurais dû te faire passer ce lien : http://www.fullerton.edu/math/mathews/a2001/Animations/RootFinding/NewtonMethod/NewtonMethod.html Avec les animations, tu devrais bien comprendre ce qui se passe, au moins pour les fonctions d'une variable. Je vois pas l'intérêt de coder Newton-Raphson sous Matlab alors qu'elle l'ai déjà. Et les fonctions de recherche de zéros de Matlab sont certainement mieux implémentées et plus efficaces qu'une méthode de Newton-Raphson codée soi-même (choix de la méthode dans les fonctions de recherche de zéros, options diverses, etc, etc...). Mais ça prend pas longtemps à en coder une si ça te fait plaisir... Faudra peut-être enlever le truc de l'erreur sur la solution elle-même, y a sûrement des cas pathologiques pour lesquels l'update est petit alors qu'on est encore loin de la solution. function [xsol errx errf]=NR(f,x0,errxmax,errfmax,nmax,fporh) %% f is a function handle %% x0 is the initial guess %% errxmax is the desired maximum estimated error on the root of f %% Put errxmax=Inf to search the root of f without caring about it %% errfmax is the desired maximum error on f(xsol) %% Put errfmax=Inf to search the root of f without caring about it %% nmax is the maximum number of iterations %% Put nmax=Inf to search the root %% without carrying about the maximum %% number of iterations %% Dangerous in pathologic cases, the while loop %% can never end %% fporh is %% either a function handle : the derivative of f %% or real number : the step for the numerical computation %% of the derivative of f xsol=x0; errf=abs(f(x0)); n=0; %% if isa(fporh,'numeric') h=fporh; errx=abs(-2*h*f(x0)/(f(x0+h)-f(x0-h))); while (n=errxmax || errf>=errfmax ... || (errxmax==Inf && errfmax==Inf))) updatexsol=-2*h*f(xsol)/(f(xsol+h)-f(xsol-h)); xsol=xsol+updatexsol; n=n+1; errx=abs(-2*h*f(xsol)/(f(xsol+h)-f(xsol-h))); errf=abs(f(xsol)); end elseif isa(fporh,'function_handle') fp=fporh; errx=abs(-f(xsol)/fp(xsol)); while (n=errxmax || errf>=errfmax ... || (errxmax==Inf && errfmax==Inf))) updatexsol=-f(xsol)/fp(xsol); xsol=xsol+updatexsol; n=n+1; errx=abs(-f(xsol)/fp(xsol)); errf=abs(f(xsol)); end else disp('The last argument has to be') disp('either a function handle or a real number') error('Error : try again...') end Ciao

Sacabouffe · Answer

Salut Au cas où t'aurais la Symbolic Math Toolbox aussi, tu peux te passer de donner la dérivée ou le pas comme argument pour NR : function [xsol errx errf]=NR(f,x0,errxmax,errfmax,nmax) %% f is a symbolic expression %% x0 is the initial guess %% errxmax is the desired maximum estimated error on the root of f %% Put errxmax=Inf to search the root of f without caring about it %% errfmax is the desired maximum error on f(xsol) %% Put errfmax=Inf to search the root of f without caring about it %% nmax is the maximum number of iterations %% Put nmax=Inf to search the root %% without carrying about the maximum %% number of iterations %% Dangerous in pathologic cases, the while loop %% can never end x=findsym(f); fp=diff(f,x); %% xsol=x0; errf=abs(subs(f,x,x0)); n=0; %% errx=abs(-subs(f,x,xsol)/subs(fp,x,xsol)); while (n=errxmax || errf>=errfmax ... || (errxmax==Inf && errfmax==Inf))) updatexsol=-subs(f,x,xsol)/subs(fp,x,xsol); xsol=xsol+updatexsol; n=n+1; errx=abs(-subs(f,x,xsol)/subs(fp,x,xsol)); errf=abs(subs(f,x,xsol)); end Ou alors tu rajoutes ce petit morceau au code du message 12. Dans ce cas, ça donne un truc de ce style : function [xsol errx errf]=NR(f,x0,errxmax,errfmax,nmax,fporh) %% f is a function handle or a symbolic expression %% x0 is the initial guess %% errxmax is the desired maximum estimated error on the root of f %% Put errxmax=Inf to search the root of f without caring about it %% errfmax is the desired maximum error on f(xsol) %% Put errfmax=Inf to search the root of f without caring about it %% nmax is the maximum number of iterations %% Put nmax=Inf to search the root %% without carrying about the maximum %% number of iterations %% Dangerous in pathologic cases, the while loop %% can never end %% fporh has to be given if f isn't a symbolic expression %% In this case, it is : %% either a function handle : the derivative of f %% or real number : the step for the numerical computation %% of the derivative of f n=0; %% if isa(f,'function_handle') xsol=x0; errf=abs(f(x0)); if isa(fporh,'numeric') h=fporh; errx=abs(-2*h*f(x0)/(f(x0+h)-f(x0-h))); while (n=errxmax || errf>=errfmax ... || (errxmax==Inf && errfmax==Inf))) updatexsol=-2*h*f(xsol)/(f(xsol+h)-f(xsol-h)); xsol=xsol+updatexsol; n=n+1; errx=abs(-2*h*f(xsol)/(f(xsol+h)-f(xsol-h))); errf=abs(f(xsol)); end elseif isa(fporh,'function_handle') fp=fporh; errx=abs(-f(xsol)/fp(xsol)); while (n=errxmax || errf>=errfmax ... || (errxmax==Inf && errfmax==Inf))) updatexsol=-f(xsol)/fp(xsol); xsol=xsol+updatexsol; n=n+1; errx=abs(-f(xsol)/fp(xsol)); errf=abs(f(xsol)); end else disp('The last argument has to be') disp('either a function handle or a real number') error('Error : try again...') end elseif isa(f,'sym') x=findsym(f); fp=diff(f,x); xsol=x0; errf=abs(subs(f,x,x0)); errx=abs(-subs(f,x,xsol)/subs(fp,x,xsol)); while (n=errxmax || errf>=errfmax ... || (errxmax==Inf && errfmax==Inf))) updatexsol=-subs(f,x,xsol)/subs(fp,x,xsol); xsol=xsol+updatexsol; n=n+1; errx=abs(-subs(f,x,xsol)/subs(fp,x,xsol)); errf=abs(subs(f,x,xsol)); end else disp('The first argument has to be') disp('either a function handle or symbolic expression') error('Error : try again...') end Bon après-midi

atiffa · Answer

Merci sacabouffe  :-)

J'ai beau essayé de comprendre mais je n'y parvient tout de même pas, qu'est ce qu'il est compliqué cette méthode (j'ai compris avec l'animation de ce lien : http://www.fullerton.edu/math/ )
mais j'ai presque rien compris dans le code que tu ma envoyer

Mais j'ai compris que tu es une personne bien, merci pour ton aide c'est très gentil 

Bon après-midi à toi aussi


Cordialement
atiffa

Sacabouffe · Answer

Y a pas de quoi ;-)

Je veux bien essayer de résumer pour accompagner l'explication graphique, au cas où. Pour une fonction d'une variable ça se passe comme ça (c'est le même principe pour des fonctions de plusieurs variables mais c'est peut-être plus dur à se représenter, dans ce cas, ça marche avec le gradient) :
On prend un initial guess x0
Ensuite le processus itératif est le suivant
1) On calcule l'équation de la tangente à la courbe au point d'abscisse x0, c'est-à-dire au point de coordonnées (x0,(f(x0))
2) Le réel x1 de l'itération suivante dans la méthode de Newton-Raphson est l'abscisse du point d'intersection de cette tangente avec l'axe des abscisses
3) Tu réitères le procédé autant de fois que tu veux ou jusqu'à atteindre un certaine précision pour f(x0) ou encore jusqu'à ce que la différence obtenue entre deux itérations vaille pas la peine d'aller plus loin.

C'est ce que fait le code.

T'as trois possibilités pour passer la fonction en argument.
1) Tu donnes la fonction et la dérivée sous forme de function handle
2) Tu donnes la fonction sous forme de function handle et un pas pour calculer numériquement la dérivée à l'aide de la formule de Taylor au premier ordre
3) Tu donnes seulement la fonction sous forme d'expression symbolique et la dérivée est calculée symboliquement par Matlab

Bonne nuit

atiffa · Answer

Merci sacabouffe,

J'avoue que c'est un peu compliquer pour ma petite tête de comprendre tous ca...
Merci encore je vais essayer de comprendre et j'espère arriver a la solution


Merci et bonne nuit


Cordialement
atiffa

atiffa · Answer

Bonjour sacabouffe :-)

Justement ce qu'il me prenne la tête dans mon code c'est la partie physique (Le courant ionique et le potentiel de la membrane), parce que le flux des ions à travers la membrane cellulaire est entraîné par des gradients de concentration et par le potentiel électrique ( on représente toujours la membrane cellulaire comme un circuit électrique)

Tu vois qu'on peut pas séparer les deux domaines ;-)))))

Bonne soirée et bonne WE


Amicalement
atiffa

[Matlab] Méthode de Newton-Raphson

16 réponses

Discussions similaires