1

MATLAB で Newton-Raphson 法を使用して、3 変数の 3 つの非線形システムを解こうとしています。3 つの非線形方程式を次に示します。

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

ニュートン・ラフソン法を使用して、、、cおよびsの値を見つける必要があります。q

これが私がこれまでに持っているものです：

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

誰かが私のコードをチェックしてください。いくつかのエラーがあるため、動作していません。前もって感謝します。

4

1 に答える 1

1

c、s、およびq( x0) の推測を反復間で更新していません。次のことを試してください。

%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

上記のように変更x0しxold、ループの繰り返しごとに更新されるようにしました。基本的に、この方法を使用して、推測を公称値に近づけようとしています。このアルゴリズムに関連する基本的な概念については、この Web サイトを参照してください。

上記の変更により、コードは 6 回の反復後に次のように収束します。

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
于 2014-06-17T20:02:13.240 に答える