Newton Raphson to solve a system of non-linear equations

47 views (last 30 days)
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
  1 Comment
Matt J
Matt J on 2 Dec 2023
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);

Sign in to comment.

Answers (1)

Matt J
Matt J on 2 Dec 2023
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);

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!