applying while loop for solving simultaneous equations

1 次查看（过去 30 天）
Akshay Pratap Singh
I am trying to include a condition that 'tan(alphac)<(1/lam)' in the following code. Here, ''alphac'' is dependent of ''x''. and the value of x is determined at last. I want to write a code such that the value of x is determined by applying the condition.
code:
dbstop if error
clear all
clc
format longEng
syms x y lam
a=[4;0.55];
% The Newton-Raphson iterations starts here
LAM=linspace(0,10,11);
h=4;
q=20;
gma=18.4; nq=2*q/(gma*(h+x));
delta=26;
phi=39;
A=lam*nq/(1+nq);
kv=0;
kh=0;
da1=delta*(pi/180); da2=-delta*(pi/180); pha1=phi*(pi/180); pha2=phi*(pi/180);
dp1=delta*(pi/180); dp2=delta*(pi/180); php1=phi*(pi/180); php2=phi*(pi/180);
psi=atan(kh/(1-kv));
m=pha1+da1;
b=pha1-psi;
c=psi+da1;
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-pha1)+(kh/(1-kv)))/(tan(alphac)*(cos(da1)+sin(da1)*tan(alphac-pha1)));
r=1-lam*tan(alphac);
kq=r*kg;
A2=0;
alphac2=atan((sin(m)*sin(b)+(sin(m)^2*sin(b)^2+sin(m)*cos(m)*sin(b)*cos(b)+A2*cos(c)*cos(m)*sin(b))^0.5)/(A2*cos(c)+sin(m)*cos(b)))
kg2=(tan(alphac2-pha1)+(kh/(1-kv)))/(tan(alphac2)*(cos(da1)+sin(da1)*tan(alphac2-pha1)));
pg=0.5*gma*(1-kv)*kg2*(h+x)^2;
A=1;
alphac1=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)))
kg1=(tan(alphac1-pha1)+(kh/(1-kv)))/(tan(alphac1)*(cos(da1)+sin(da1)*tan(alphac1-pha1)));
r1=1-lam*tan(alphac1);
kq1=r1*kg1;
pq=(1-kv)*(q*kq*(h+x)+0.5*q*(kq1-kq)*(h+x));
va2=asin(sin(da2)/sin(pha2))-asin(sin(psi)/sin(pha2))-da2-psi;
ka2=(1/cos(psi))*(cos(da2)*((cos(da2)-sqrt(sin(pha2)^2-sin(da2)^2)))/(cos(psi)+sqrt(sin(pha2)^2-sin(psi)^2)))*exp(-va2*tan(pha2));
vp1=asin(sin(dp1)/sin(php1))+asin(-sin(psi)/sin(php1))+dp1+psi;
kp1=(1/cos(psi))*(cos(dp1)*((cos(dp1)+sqrt(sin(php1)^2-sin(dp1)^2)))/(cos(psi)-sqrt(sin(php1)^2-sin(psi)^2)))*exp(vp1*tan(php1));
vp2=asin(sin(dp2)/sin(php2))+asin(-sin(psi)/sin(php2))+dp2+psi;
kp2=(1/cos(psi))*(cos(dp2)*((cos(dp2)+sqrt(sin(php2)^2-sin(dp2)^2)))/(cos(psi)-sqrt(sin(php2)^2-sin(psi)^2)))*exp(vp2*tan(php2));
sinda1=sin(da1); sindp1=sin(dp1); sinda2=-sin(da2); sindp2=sin(dp2);
cosda1=cos(da1); cosdp1=cos(dp1); cosda2=cos(da2); cosdp2=cos(dp2);
pp1=kp1*gma*0.5*(x^2);
pa2=ka2*gma*(x*y+0.5*(y^2)); pp2=kp2*gma*(y*(h+x)+(0.5*(y^2)));
zp1=x/3;
zp2=((0.5*(h+x)*(y^2))+((y^3)/3))/(((h+x)*y)+(0.5*(y^2)));
za2=((0.5*x*(y^2))+((y^3)/3))/((x*y)+(0.5*(y^2)));
e2=(pp1*cosdp1)+(pa2*cosda2)-(pg*cosda1)-(pp2*cosdp2)-pq;
e3=(pp1*cosdp1*zp1)+(pp2*cosdp2*zp2)-(pg*cosda1*((h+x)/3))-(pa2*cosda2*za2)-pq*(1/3)*(h+x)*((kq+2*kq1)/(kq+kq1));
g=[e2; e3];
J=jacobian([e2, e3], [x, y]);
A=zeros(2,numel(LAM));
for i=1:numel(LAM)
del=1;
indx=0;
lam=0;
while del>1e-6 && tan(alphac)<(1/lam)
gnum = vpa(subs(g,[x,y,lam],[a(1),a(2),LAM(i)]));
Jnum = vpa(subs(J,[x,y,lam],[a(1),a(2),LAM(i)]));
delx = -Jnum\gnum;
a = a + delx;
del = max(abs(gnum));
indx = indx + 1;
end
Z(:,i)=double(a)
end

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

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by