Méthode Newton-Raphson pour un système non linéaire de 3 variables dans Matlab

Question

J'essaie de résoudre 3 systèmes non linéaires de 3 variables en utilisant la méthode Newton-Raphson dans MATLAB. Voici les 3 équations non-linéaires:

c * (alpha*I + k_f + k_d + k_n * s + k_p*(1-q))-I *alpha  = 0 

s * (lambda_b * c* P_C + lambda_r *(1-q))- lambda_b* c * P_C = 0 

q * (gamma + c * k_p *(P_C/P_Q))- (c * k_p * (P_C/P_Q)) = 0 

Je dois trouver les valeurs de c, s et q en utilisant la méthode de Newton-Raphson.

Voici ce que j'ai jusqu'à présent:

format long 
clear; 

%values of parameters 
I=1200; 
k_f= 6.7*10.^7; 
k_d= 6.03*10.^8; 
k_n=2.92*10.^9; 
k_p=4.94*10.^9; 
lambda_b= 0.0087; 
lambda_r =835; 
gamma =2.74; 
alpha =1.14437*10.^-3; 
P_C= 3 * 10.^(11); 
P_Q= 2.87 * 10.^(10); 

tol = 10.^-4; %tol is a converge tolerance 

%initial guess or values 
c=1; 
s=0.015; 
q=0.98; 
x0= [c;s;q]; 

iter= 0; %iterations 
xnew =[100;100;100]; 
while norm(xnew -x0) > tol 
    iter= iter + 1; 
%Defining the functions for c,s and q. 
f = c * (alpha*I + k_f + k_d + k_n * s + k_p*(1-q))-I *alpha; 
g = s * (lambda_b * c* P_C + lambda_r *(1-q))- lambda_b* c * P_C; 
h = q * (gamma + c * k_p *(P_C/P_Q))- (c * k_p * (P_C/P_Q)); 

%Partial derivatives in terms of c,s and q. 
dfdc = alpha*I + k_f + k_d + k_n * s + k_p*(1-q); 
dfds = k_n *c ; 
dfdq = - k_p *c; 

dgdc = lambda_b * P_C *(s-1); 
dgds = lambda_b * c* P_C + lambda_r *(1-q); 
dgdq = - lambda_r * s; 

dhdc = k_p *(P_C/P_Q)*(q-1); 
dhds = 0; 
dhdq = gamma + c * k_p *(P_C/P_Q); 

%Jacobian matrix 
J = [dfdc dfds dfdq; dgdc dgds dgdq; dhdc dhds dhdq];  
% Applying the Newton-Raphson method 
xnew = x0 - J\[f;g;h]; 
disp(sprintf('iter=%6.15f, c=%6.15f, s=%6.15f, q=%6.15f', iter,xnew)); 
end 

quelqu'un peut-il s'il vous plaît vérifier mon code, il y a quelques erreurs donc, sa ne fonctionne pas. Merci d'avance.

Answer 1

Vous ne mettez pas à jour votre estimation pour c, s et q (x0) entre les itérations. Effectuez les opérations suivantes:

%initial guess or values 
c=1; 
s=0.015; 
q=0.98; 

xnew =[c;s;q]; 
xold = zeros(size(xnew)); 
while norm(xnew - xold) > tol 
    iter= iter + 1; 
    xold = xnew; 

    % update c, s, and q 
    c = xold(1); 
    s = xold(2); 
    q = xold(3); 

    %Defining the functions for c,s and q. 
    f = c * (alpha*I + k_f + k_d + k_n * s + k_p*(1-q))-I *alpha; 
    g = s * (lambda_b * c* P_C + lambda_r *(1-q))- lambda_b* c * P_C; 
    h = q * (gamma + c * k_p *(P_C/P_Q))- (c * k_p * (P_C/P_Q)); 

    %Partial derivatives in terms of c,s and q. 
    dfdc = alpha*I + k_f + k_d + k_n * s + k_p*(1-q); 
    dfds = k_n *c ; 
    dfdq = - k_p *c; 

    dgdc = lambda_b * P_C *(s-1); 
    dgds = lambda_b * c* P_C + lambda_r *(1-q); 
    dgdq = - lambda_r * s; 

    dhdc = k_p *(P_C/P_Q)*(q-1); 
    dhds = 0; 
    dhdq = gamma + c * k_p *(P_C/P_Q); 

    %Jacobian matrix 
    J = [dfdc dfds dfdq; dgdc dgds dgdq; dhdc dhds dhdq];  
    % Applying the Newton-Raphson method 
    xnew = xold - J\[f;g;h]; 
    disp(sprintf('iter=%6.15f, c=%6.15f, s=%6.15f, q=%6.15f', iter,xnew)); 
end 

j'ai changé x0 être xold dans ce qui précède, et fait en sorte qu'il est mis à jour à chaque itération de la boucle. Fondamentalement, vous essayez de "marcher" vos estimations plus près des valeurs nominales en utilisant cette méthode. Découvrez this website pour une bonne discussion sur les concepts de base impliqués dans cet algorithme.

Avec les modifications ci-dessus, la fin de votre code de configuration convergeant après six itérations avec:

iter=1.000000000000000, c=0.000000000389029, s=0.015000000287216, q=0.979999999955780 
iter=2.000000000000000, c=0.000000001356998, s=0.158028331191731, q=0.923765241962920 
iter=3.000000000000000, c=0.000000001181617, s=0.104156261426515, q=0.952886937302707 
iter=4.000000000000000, c=0.000000001216663, s=0.085849634576983, q=0.958360887077671 
iter=5.000000000000000, c=0.000000001224388, s=0.084367460093642, q=0.958463129596494 
iter=6.000000000000000, c=0.000000001224423, s=0.084367992582976, q=0.958463625488450