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

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

Answer 1

[Torsten]

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.

[Yamini]

how can I set constraint in the code . can you please elaborate or can give an example code for it.

[Torsten]

If you look at your equations, you will notice that [0,0,0] is a solution. And it seems that the solver always converges to this point.

z0 = [3; 1.646; 33.754];  % Initial guess for [x; y; P]
lb = []; % set lower bounds in [x,y,P]
ub = []; % set upper bounds on [x,y,P]
options = optimset('TolX',1e-7,'TolFun',1e-7);
sol = lsqnonlin(@fun,z0,lb,ub,[],[],[],[],[],options)

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

sol = 3x1

0.0017 -0.0002 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fun(sol)

ans = 1x3

1.0e-04 * 0.7326 0.1049 0.0025

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function res = fun(z)

syms x y P y1 y2

% Define constants
phi = (pi / 180) * 30;
delta = (2 * phi) / 3;
gma = 14.2;

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;
VF = -v1 + v2 + v3 - v4;
MF = m2 + m4 - m1 - m3 + P * (h + x);
res = double(subs([HF,VF,MF],[x,y,P],[z(1),z(2),z(3)]));
end

[Yamini]

Can you please write the full code once by making changes in my initial code given in question because I am not getting how to set constraint in the code. Also I think that P_vals = linspace(0, 33.754, 10); is not making any sense because P is unknown variable instead we can do it for surcharge,q as Q= linspace(l1.l2,8);. so please help me in writing the code.

[Torsten]

I copied the part of your code in which you compute HF, VF and MF into the function fun and called "lsqnonlin" on this problem. So the code is complete.

The bounds on the parameters are to be set in the arrays lb and ub above.

[Yamini]

can you please write the full code so that it works directly on my matlab and I can see is this process is giving right results.

[Torsten]

This is the full code which works in MATLAB R2024a. You only need to insert the bounds for x, y and P in lb and ub.

If the code does not work in your MATLAB version: what is your MATLAB version and which error message do you get ? Do you have a license for the Optimization Toolbox ?

z0 = [3; 1.646; 33.754];  % Initial guess for [x; y; P]
lb = []; % set lower bounds in [x,y,P]
ub = []; % set upper bounds on [x,y,P]
options = optimset('TolX',1e-7,'TolFun',1e-7);
sol = lsqnonlin(@fun,z0,lb,ub,[],[],[],[],[],options)
fun(sol)

function res = fun(z)

syms x y P y1 y2

% Define constants
phi = (pi / 180) * 30;
delta = (2 * phi) / 3;
gma = 14.2;

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;
VF = -v1 + v2 + v3 - v4;
MF = m2 + m4 - m1 - m3 + P * (h + x);
res = double(subs([HF,VF,MF],[x,y,P],[z(1),z(2),z(3)]));
end

[Yamini]

I am using MATLAB online version. I am getting the result but they are not correct instead the value of the variable should be as initial guess or with a variation within the range of 10%-15% from the initial guess. I think I have license for Optimization Toolbox. So can you make the required changes so that it work efficiently for my version also.

[Torsten]

I am getting the result but they are not correct instead the value of the variable should be as initial guess or with a variation within the range of 10%-15% from the initial guess.

Then set

z0 = [3; 1.646; 33.754];  % Initial guess for [x; y; P]
lb = z0*0.85; % set lower bounds in [x,y,P]
ub = z0*1.15; % set upper bounds on [x,y,P]

and see what you get in the 15% region around your initial guess values. But I think that the solver will converge to lb because it will head towards [0;0;0] which solves your problem.

So can you make the required changes so that it work efficiently for my version also.

You still didn't tell us what problems you have with your MATLAB version.

[Yamini]

The code seems like we are just taking out the result as exactly the 15% +- initial guess. But i dont want it. The answer should be coming automatically within the range from the initial guess to +-(10%-15%). i think you should once check the initial code given by me in the question and should suggest the changes I should make but not by changing the entire method. I had got the results from newton rapson iteration method using jacobian matrix for solving equation with 2 variables but I am finding difficulty for 3 variables. I think something need to be corrected at the last part after the starting of newton rapson iteration method.

and there is no problem with my matlab version its just that I dont have version which yu were asking me about.

[Torsten]

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.

[Yamini]

okay thankyou for help