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

Question

Yamini 2024-9-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2150329-give-correction-to-my-code-for-finding-value-of-three-variable-by-newton-rapson-iteration

评论： Yamini 2024-9-11

在 MATLAB Online 中打开

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

% 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

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-9-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2150329-give-correction-to-my-code-for-finding-value-of-three-variable-by-newton-rapson-iteration#answer_1511389

编辑：Torsten 2024-9-5

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 个评论
显示 9更早的评论隐藏 9更早的评论

Torsten 2024-9-5

编辑：Torsten 2024-9-5

在 MATLAB Online 中打开

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

Torsten 2024-9-9

编辑：Torsten 2024-9-9

在 MATLAB Online 中打开

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

Torsten 2024-9-11

编辑：Torsten 2024-9-11

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 2024-9-11

okay thankyou for help

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

11 个评论
显示 9更早的评论隐藏 9更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

11 个评论 显示 9更早的评论隐藏 9更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

11 个评论
显示 9更早的评论隐藏 9更早的评论