solution by newton raphson method for nonlinear equations

2 views (last 30 days)
I am writing a code for solving two non linear simultaneous equations using newton raphson method. I am not able to link the g and J for different variables with newton raphson method. As I am new to matlab. Please help and thank in advance.
dbstop if error
clear all
clc
format longEng
syms x y h
phi=(pi/180)*39;
delta=(2*phi)/3;
gma=18.4;
a=[1.5;0.632];
% kh=0.3;
H=linspace(0.5,4,8);
% h=4;
lam=0;
q=50;
nq=2*q/(gma*(h+x));
A=lam*nq/(1+nq);
kh=0;
kv=0;
psi=atan(kh/(1-kv));
beta=1.2;
alfa=1.2;
R1=-1;
R3=-1;
delm1=0.5*(1-R1)*delta;
delm3=-0.5*(1-R3)*delta;
m=phi+delm1;
b=phi-psi;
c=psi+delm1;
alphac=atan((sin(m)*sin(b)+(sin(m)^2*sin(b)^2+sin(m)*cos(m)*sin(b)*cos(b)+A*cos(c)*cos(m)*sin(b))^0.5)/(A*cos(c)+sin(m)*cos(b)));
kg=(tan(alphac-phi)+(kh/(1-kv)))/(tan(alphac)*(cos(delm1)+sin(delm1)*tan(alphac-phi)));
r=1-lam*tan(alphac);
kq=r*kg;
pg=0.5*gma*(1-kv)*kg*(h+x)^2*cos(delm1);
pq=(1-kv)*q*kq*(h+x);
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);
y1=0:0.01:1;
y2=0:0.01:1;
mmm=size(y1);
mm=mmm(1,2);
for i=1:mm
if (y1(i)>=1-(1/beta)&& y1(i)<=1)
R2(i)=3*(beta*(1-y1(i))).^0.5;
else
R2(i)=3;
end
if (R2(i)>=0 && R2(i)<1)
delm2(i)=0.5*(1-R2(i)).*delta;
k2(i)=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R2(i))+cos(psi)*cos(delm2(i)+psi)...
*(1-R2(i))*(1+sqrt((sin(phi+delm2(i))*sin(phi-psi))/cos(delm2(i)+psi)))^2);
else
delm2(i)=0.5*(R2(i)-1)*delta;
k2(i)=1+0.5*(R2(i)-1).*((cos(phi-psi).^2./(cos(psi).*(cos(delm2(i)+psi).*...
(-sqrt((sin(phi+delm2(i)).*sin(phi-psi))./(cos(delm2(i)+psi)))+1).^2)))-1);
end
if (y2(i)>=0 && y2(i)<=(1/alfa))
R4(i)=3*(alfa*(y2(i)))^0.5;
else
R4(i)=3;
end
if (R4(i)>=0 && R4(i)<=1)
delm4(i)=0.5*(1-R4(i)).*delta;
k4(i)=(2*cos(phi-psi)^2)/(cos(phi-psi)^2*(1+R4(i))+cos(psi)*cos(delm4(i)+psi)...
*(1-R4(i))*(1+sqrt((sin(phi+delm4(i))*sin(phi-psi))/cos(delm4(i)+psi)))^2);
else
delm4(i)=0.5*(R4(i)-1)*delta;
k4(i)=1+0.5*(R4(i)-1)*((cos(phi-psi)^2/(cos(psi)*(cos(delm4(i)+psi)*(-sqrt((sin(phi+delm4(i))...
*sin(phi-psi))/(cos(delm4(i)+psi)))+1)^2)))-1);
end
h2(i)=gma.*x.^2.*k2(i).*y1(i).*cos(delm2(i));
h4(i)=gma.*y.^2.*k4(i).*cos(delm4(i)).*y2(i)+gma.*(x+h).*y.*k4(i).*cos(delm4(i));
H3_a(i)=k3.*y2(i).*cos(delm3);
% H3_b=matlabFunction(k3*cos(delm3));
h3=gma.*y.^2.*H3_a+gma.*x.*y.*k3.*cos(delm3);%integral(H3_b,0,1);
HF=h2-h4+h3-pg-pq;
%
M2(i)=k2(i).*y1(i).*cos(delm2(i)).*(1-y1(i));
m2=gma.*x.^3.*M2;
M4_a(i)=y2(i).^2;
M4_b(i)=y2(i);
m4=gma.*y.^3.*k4(i).*cos(delm4(i)).*M4_a+gma.*(x+h).*y.^2.*k4(i).*cos(delm4(i)).*M4_b;
M3_a(i)=k3.*y2(i).^2.*cos(delm3);
M3_b(i)=k3.*y2(i).*cos(delm3);
m3=gma.*y.^3.*M3_a+gma.*x.*y.^2.*M3_b;
MF=m2+m4-m3-pg.*(h+x).*(1/3)-0.5.*pq.*(h+x);
g=[HF; MF];
J=jacobian([HF, MF], [x, y]);
end
% Z=zeros(2,numel(H));
for j=1:numel(H)
del=1;
indx=0;
while del>1e-6
gnum = vpa(subs(g,[x,y],[a(1),a(2)]));
Jnum = vpa(subs(J,[x,y],[a(1),a(2)]));
delx = -Jnum.\gnum;
a = a + delx;
del = max(abs(gnum));
indx = indx + 1;
end
Z(:,j)=double(a)
end
'NEWTON-RAPHSON SOLUTION CONVERGES IN ITERATIONS',indx,
'FINAL VALUES OF a ARE';a,

Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!