Newton Raphson to solve a system of non-linear equations

I am trying to solve a system of non-linear equations using Newton-Raphson using the code I listed at the bottom.

I tried matlab's fsolve and that gave me the solution

0.0136 1.1846

My program seem to give the wrong solution and give me a warning for an ill conditioned matrix if I decrease the error/tolerance.

I don't understand what I am doing wrong. I have already checked the functions I am using, and tried to rely on by-hand derivation of the Jacobian rather than the built in matlab function, but something seems to be off.

Any help is greatly appreciated.

I am solving for d and D. I set eq(1) = eq(2) and plugged in the reciprocal of U_f to get the first equation

the third equation is function (3) - delta P = 29.6

clear all;
syms d D


%Constants
R_fi = 1.76e-4;     R_fo = 1.76e-4; 
h_s = 356;          h_t = 356; 
k_w = 60;           c = 0.389;
S_t = 0.016;        K_1 = 0.249;
n_1 = 2.207;        Q = 801368;

delta_T_m = 29.6;   D_s = 1.219; 
delta_P = 49080 ;   A_ex_condition = 64.15; 

%__________________________________________________________________________

U_f_reciprocal = (((((d/(D*h_t)) + (d*R_fi/D) + (d*log(d/D)/(2*k_w)) + R_fo + (1/h_s)))));

f_1 = (pi * K_1 * (D/d))-(Q * (1/delta_T_m) * U_f_reciprocal);
f_2 = (((D^2) * ((S_t-d)^2))/c)-(1/delta_P);

F_express = [f_1;f_2];
Df_express = jacobian([f_1,f_2],[d D]);

Df = matlabFunction(Df_express);
F = matlabFunction(F_express);
%__________________________________________________________________________
x = [0.015;0.8]; % starting guess

err=1;
tol=10e-2;


while err>tol
    Dx = -1 * Df(x(1),x(2))\(F(x(1),x(2))); % solve for increment
    x = x + Dx; % add on to get new guess
    F(x(1),x(2)); 
    err = abs(x-Dx);
end

[Matt J]

Did you mean for the err variable in your while loop to be a scalar? If so, it should be something like,

err = norm(x-Dx,1);

Answer 1

[Matt J]

From the surface plot, the equations don't appear to have any solution, or at least no where near the vicinity you are searching. Also, Newton-Raphson is not really well suited to this problem because it doesn't allow you to confine the search to positive d and D, which appears to be necessary in order for the log expressions in your equations to make sense.

syms d D

%Constants
R_fi = 1.76e-4;     R_fo = 1.76e-4; 
h_s = 356;          h_t = 356; 
k_w = 60;           c = 0.389;
S_t = 0.016;        K_1 = 0.249;
n_1 = 2.207;        Q = 801368;

delta_T_m = 29.6;   D_s = 1.219; 
delta_P = 49080 ;   A_ex_condition = 64.15; 

%__________________________________________________________________________

U_f_reciprocal = (((((d/(D*h_t)) + (d*R_fi/D) + (d*log(d/D)/(2*k_w)) + R_fo + (1/h_s)))));

f_1 = (pi * K_1 * (D/d))-(Q * (1/delta_T_m) * U_f_reciprocal);
f_2 = (((D^2) * ((S_t-d)^2))/c)-(1/delta_P);

F_express = [f_1;f_2];
Df_express = jacobian([f_1,f_2],[d D]);

Df = matlabFunction(Df_express);
F = matlabFunction(F_express);

fsurf(@(d,D) norm(F(d,D))  ,[0.01,0.2,.1,5]); xlabel d; ylabel D

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

view(35,20);