Clear Filters
Clear Filters

give correction to my code for finding value of three variable by newton rapson iteration.

9 views (last 30 days)
my matlab code:
clear all
clc
format longEng
syms x y P y1 y2
% Define constants
phi = (pi / 180) * 30;
delta = (2 * phi) / 3;
gma = 14.2;
t = [3; 1.646]; % Initial guess for x and y
P_vals = linspace(0, 33.754, 10);
h = 5;
q = 0;
a = 0;
kh = 0;
kv = 0;
psi = atan(kh / (1 - kv));
beta = 6;
alfa = 0.7;
R1 = -1;
R3 = -1;
delm1 = 0.5 * (1 - R1) * delta;
delm3 = -0.5 * (1 - R3) * delta;
l = phi + delta;
m = phi - psi;
n = psi + delta;
mq = 2 * q / (gma * x);
A = a * mq / (x * (1 + mq));
thetacm = atan((sin(l) * sin(m) + sqrt(sin(l)^2 * sin(m)^2 + A * cos(l) * sin(m) * cos(n) + cos(m) * sin(l) * sin(m) * cos(l))) / (cos(m) * sin(l) - A * cos(n)));
B = (x / tan(thetacm)) - a;
kg = (tan(thetacm - phi) + tan(psi)) / (tan(thetacm) * (cos(2 * phi / 3) + sin(2 * phi / 3) * tan(thetacm - phi)));
kq = (1 - atan(thetacm) / x) * ((tan(thetacm - phi) + tan(psi)) / (tan(thetacm) * (cos(2 * phi / 3) + sin(2 * phi / 3) * tan(thetacm - phi))));
pg = 0.5 * gma * (1 - kv) * kg * x^2 * cos(delm1);
pq = (1 - kv) * q * kq * (B / (B + 2 * a)) * (x - a) * cos(delm1);
k4R0 = (2 * cos(phi - psi)^2) / (cos(phi - psi)^2 + cos(psi) * cos(delta / 2 + psi) * (1 + sqrt(sin(phi + delta / 2) * sin(phi - psi) / cos(delta / 2 + psi)))^2);
k3 = (2 * cos(phi - psi)^2) / (cos(phi - psi)^2 * (1 + R3) + cos(psi) * cos(delm3 + psi) * (1 - R3) * (1 + sqrt(sin(phi + delm3) * sin(phi - psi) / cos(delm3 + psi)))^2);
delm23 = 0.5 * (3 - 1) * delta;
k23 = 1 + 0.5 * (3 - 1) * ((cos(phi - psi)^2 / (cos(psi) * (cos(delm23 + psi) * (-sqrt(sin(phi + delm23) * sin(phi - psi) / cos(delm23 + psi))) + 1)^2)) - 1);
R2 = 3 * (beta * (1 - y1))^0.5;
delm213 = 0.5 * (R2 - 1) * delta;
k213 = 1 + 0.5 * (R2 - 1) * ((cos(phi - psi)^2 / (cos(psi) * (cos(delm213 + psi) * (-sqrt(sin(phi + delm213) * sin(phi - psi) / cos(delm213 + psi))) + 1)^2)) - 1);
delm201 = 0.5 * (1 - R2) * delta;
k201 = (2 * cos(phi - psi)^2) / (cos(phi - psi)^2 * (1 + R2) + cos(psi) * cos(delm201 + psi) * (1 - R2) * (1 + sqrt(sin(phi + delm201) * sin(phi - psi) / cos(delm201 + psi)))^2);
delm43 = 0.5 * (3 - 1) * delta;
k43 = 1 + 0.5 * (3 - 1) * ((cos(phi - psi)^2 / (cos(psi) * (cos(delm43 + psi) * (-sqrt(sin(phi + delm43) * sin(phi - psi) / cos(delm43 + psi))) + 1)^2)) - 1);
R4 = 3 * (alfa * y2)^0.5;
delm413 = 0.5 * (R4 - 1) * delta;
k413 = 1 + 0.5 * (R4 - 1) * ((cos(phi - psi)^2 / (cos(psi) * (cos(delm413 + psi) * (-sqrt(sin(phi + delm413) * sin(phi - psi) / cos(delm413 + psi))) + 1)^2)) - 1);
delm401 = 0.5 * (1 - R4) * delta;
k401 = (2 * cos(phi - psi)^2) / (cos(phi - psi)^2 * (1 + R4) + cos(psi) * cos(delm401 + psi) * (1 - R4) * (1 + sqrt(sin(phi + delm401) * sin(phi - psi) / cos(delm401 + psi)))^2);
% Define integrands
H23 = matlabFunction(k23 * y1 * cos(delm23));
h23 = gma * x^2 * integral(H23, 0, (1 - (1 / beta)));
H213 = matlabFunction(k213 * cos(delm213) * y1);
h213 = gma * x^2 * integral(H213, (1 - (1 / beta)), (1 - (1 / (9 * beta))));
H201 = matlabFunction(k201 * y1 * cos(delm201));
h201 = gma * x^2 * integral(H201, (1 - (1 / (9 * beta))), 1);
h2 = h23 + h213 + h201;
H401 = matlabFunction(k401 * y2 * cos(delm401));
h401 = gma * y^2 * integral(H401, 0, (1 / (9 * alfa)));
H413 = matlabFunction(k413 * y2 * cos(delm413));
h413 = gma * y^2 * integral(H413, (1 / (9 * alfa)), (1 / alfa));
H43 = matlabFunction(k43 * y2 * cos(delm43));
h43 = gma * y^2 * integral(H43, (1 / alfa), 1);
hb43 = (x + (q / gma)) * y * k4R0 * cos(delta / 2);
h4 = h401 + h413 + h43 + hb43;
h3 = 0.5 * gma * k3 * y^2 * cos(delta) + gma * x * k3 * y * cos(delta);
h1 = pg + pq;
HF = -P - h1 + h2 + h3 - h4;
V23 = matlabFunction(k23 * y1 * sin(delm23));
v23 = gma * x^2 * integral(V23, 0, (1 - (1 / beta)));
V213 = matlabFunction(k213 * sin(delm213) * y1);
v213 = gma * x^2 * integral(V213, (1 - (1 / beta)), (1 - (1 / (9 * beta))));
V201 = matlabFunction(k201 * y1 * sin(delm201));
v201 = gma * x^2 * integral(V201, (1 - (1 / (9 * beta))), 1);
v2 = v23 + v213 + v201;
V401 = matlabFunction(k401 * y2 * sin(delm401));
v401 = gma * y^2 * integral(V401, 0, (1 / (9 * alfa)));
V413 = matlabFunction(k413 * y2 * sin(delm413));
v413 = gma * y^2 * integral(V413, (1 / (9 * alfa)), (1 / alfa));
V43 = matlabFunction(k43 * y2 * sin(delm43));
v43 = gma * y^2 * integral(V43, (1 / alfa), 1);
vb43 = (x + (q / gma)) * y * k4R0 * sin(delta / 2);
v4 = v401 + v413 + v43 + vb43;
v3 = 0.5 * gma * k3 * y^2 * sin(delta) + gma * x * k3 * y * sin(delta);
v1 = pg * cot(delm1) + pq * cot(delm1);
VF = -v1 + v2 + v3 - v4;
M23 = matlabFunction(k23 * y1 * (1 - y1) * cos(delta));
m23 = gma * x^3 * integral(M23, 0, (1 - (1 / beta)));
M213 = matlabFunction(k213 * cos(delm213) * y1 * (1 - y1));
m213 = gma * x^3 * integral(M213, (1 - (1 / beta)), (1 - (1 / (9 * beta))));
M201 = matlabFunction(k201 * y1 * cos(delm201) * (1 - y1));
m201 = gma * x^3 * integral(M201, (1 - (1 / (9 * beta))), 1);
m2 = m23 + m213 + m201;
M401 = matlabFunction(k401 * cos(delm401) * y2^2);
m401 = gma * y^3 * integral(M401, 0, (1 / (9 * alfa)));
M413 = matlabFunction(k413 * cos(delm413) * y2^2);
m413 = gma * y^3 * integral(M413, (1 / (9 * alfa)), (1 / alfa));
M43 = matlabFunction(k43 * cos(delta) * y2^2);
m43 = gma * y^3 * integral(M43, (1 / alfa), 1);
mb43 = (x + (q / gma)) * (y / 2)^2 * k4R0;
m4 = m401 + m413 + m43 + mb43;
m1 = pg * (x / 3) + pq * ((x - a) / 2);
m3 = 0.5 * gma * k3 * y^2 * cos(delta) * (2 * (y / 3)) + gma * x * k3 * y * cos(delta) * (y / 2);
MF = m2 + m4 - m1 - m3 + P * (h + x);
% Define equations
HF = -P - h1 + h2 + h3 - h4
HF = 
VF = -v1 + v2 + v3 - v4
VF = 
MF = m2 + m4 - m1 - m3 + P * (h + x)
MF = 
double(subs([HF, VF, MF],[x,y,P],[3, 1.646, 33.754]))
ans = 1x3
191.289189155057e+000 44.5300512610845e+000 534.518193150590e+000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
% Compute the Jacobian matrix
J = jacobian([HF, VF, MF], [x, y, P]);
% Initialize values
P_vals = linspace(0, 33.754, 10); % Example range for P values
Z = zeros(3, numel(P_vals)); % Matrix to store results
for i = 1:numel(P_vals)
del = 1;
indx = 0;
t = [3; 1.646; 33.754]; % Initial guess for [x; y; P]
while del > 1e-6
% Substitute current values into equations and Jacobian
subs_values = [t(1), t(2), t(3)]; % Current values for x, y, and P
gnum = double(subs([HF; VF; MF], [x, y, P], subs_values));
Jnum = double(subs(J, [x, y, P], subs_values));
% Ensure dimensions match
[m, n] = size(Jnum);
if m ~= n
error('Jacobian matrix is not square: %d x %d', m, n);
end
if length(gnum) ~= m
error('Dimension mismatch: gnum vector length (%d) does not match Jacobian matrix size (%d)', length(gnum), m);
end
% Compute update step
try
delx = -Jnum \ gnum;
catch ME
error('Matrix division failed: %s', ME.message);
end
% Check if delx is the correct size
if length(delx) ~= length(t)
error('Dimension mismatch: delx vector length (%d) does not match t vector length (%d)', length(delx), length(t));
end
% Update guess
t = t + delx;
del = max(abs(gnum));
indx = indx + 1;
end
Z(:, i) = double(t);
end
% Display results
fprintf('NEWTON-RAPHSON SOLUTION CONVERGES IN %d ITERATIONS\n', indx);
NEWTON-RAPHSON SOLUTION CONVERGES IN 16 ITERATIONS
fprintf('FINAL VALUES OF x, y, and P FOR P = %.2f ARE:\n', P_vals(end));
FINAL VALUES OF x, y, and P FOR P = 33.75 ARE:
disp(Z);
Columns 1 through 7 62.2358956951181e-006 62.2358956951181e-006 62.2358956951181e-006 62.2358956951181e-006 62.2358956951181e-006 62.2358956951181e-006 62.2358956951181e-006 -31.3808137998502e-006 -31.3808137998502e-006 -31.3808137998502e-006 -31.3808137998502e-006 -31.3808137998502e-006 -31.3808137998502e-006 -31.3808137998502e-006 1.52107049952637e-012 1.52107049952637e-012 1.52107049952637e-012 1.52107049952637e-012 1.52107049952637e-012 1.52107049952637e-012 1.52107049952637e-012 Columns 8 through 10 62.2358956951181e-006 62.2358956951181e-006 62.2358956951181e-006 -31.3808137998502e-006 -31.3808137998502e-006 -31.3808137998502e-006 1.52107049952637e-012 1.52107049952637e-012 1.52107049952637e-012
double(subs([HF, VF, MF],[x,y,P],Z(:,end).'))
ans = 1x3
1.0e+00 * 75.2337027918753e-009 10.1922481593136e-009 9.50654654042714e-012
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
this code is giving answer but the answer is not cpming as we require. We want the value of the variable coming to be as initial guess or with a variation within the range of 10%-15% from the initial guess. But this code is giving negative values and far way value from initial guess. So correct the above code so that we get the desired results. Also all the above equations are written correct I think the the part below newton rapson iteration is having some problem

Answers (1)

Torsten
Torsten on 5 Sep 2024
Edited: Torsten on 5 Sep 2024
The residual for your initial guess is [191,44,534], the residual for the last column of Z is [75e-9,10e-9,9e-12]. Thus obviously, your Newton method works, but converges to a non-physical solution. The only thing you can do is to change the initial guess or set constraints on the solution (e.g. by using "lsqnonlin" instead of your own Newton code). Further, I'd check the implemented equations for correctness.
  11 Comments
Torsten
Torsten on 11 Sep 2024
Edited: Torsten on 11 Sep 2024
The answer should be coming automatically within the range from the initial guess to +-(10%-15%).
The answer for your equations is [0;0;0], and this answer is coming with your code and with "lsqnonlin" as well. Since I don't know your model, I can't find where you possibly made a mistake or how to choose better initial guesses for the parameters to make the solver converge to the solution you have in mind. So I can't help any further.

Sign in to comment.

Products

Community Treasure Hunt

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

Start Hunting!