# plotting a graph of multivariable function

3 views (last 30 days)
Akshay Pratap Singh on 4 Dec 2019
Commented: Raj on 5 Dec 2019
I formulated a matlab code for obtaining a curve and it is giving results but inside a loop for y1 and y2. The code is taking only last values, i.e., y1 and y2=1. however, I want to include the values of y1 and y2 from 0 to 1. How can I include it.
The code is written as:
clear all
clc
format longEng
syms z y1 y2
phi=(pi/180)*39;
delta=(pi/180)*26;
gma=18.4;
h=4;
q=20;
a=0.1;
b=0.4;
h1=1.021121318;
h2=0.744439152;
L=h+h1+h2;
beta=5;
alfa=1;
% z1=z-h;
% z2=z-h-h1;
R1=-1;
R3=-1,%-(alfa*((z-h-h1)/h2))^0.5;
% R2=3*(beta*(1-((z-h)/(h1))))^0.5;
% R4=3*(alfa*((z-h-h1)/h2))^0.5;
kh=0.0;
kv=0;
psi=atan(kh/(1-kv));
da1=delta; pha1=phi;
l=pha1+da1;
m=pha1-psi;
n=psi+da1;
B=(2*q*b*(b/(2*a+b)))/(gma*(h+h1)*(h+h1));
alphacm=atan((sin(l)*sin(m)+(((sin(l)^2)*(sin(m)^2))-B*cos(n)*sin(m)*cos(l)+cos(m)*sin(l)*sin(m)*cos(l))^0.5)/((cos(m)*sin(l))-B*cos(n)))
kg1=((tan(alphacm-phi)+tan(psi))*cos(alphacm-phi))/(tan(alphacm)*(cos(alphacm-phi-delta)))
r1=(b/(h+h1))*(b/(2*a+b))*(tan(alphacm))
kq1=r1*kg1
% pgx=gma*(1-kv)*kg1*x
pqx=(1-kv)*q*kq1
%
k4R0=(2*cos(phi-psi)^2)/(cos(phi-psi)^2+cos(psi)*cos((delta/2)+psi)...
*(1+sqrt((sin(phi+(delta))*sin(phi-psi))/cos((delta/2)+psi)))^2)
%
delm3=-0.5*(1-R3)*delta;
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);
i=0;
for z=0:0.1:L
i=i+1;
if(z<=h)
p1(i)=kg1*gma*z*cos(delta)+kq1*q;
p2(i)=0;
elseif z<=(h+h1)
p1(i)=kg1*gma*z*cos(delta)+kq1*q;
for y1=0:0.01:1
if (y1>=0 && y1<=(1-(1/beta)))
p2(i)=-gma*(z-h)*k23*cos(delm23);
elseif (y1>=(1-(1/beta)) && y1<=(1-(1/(9*beta))))
p2(i)=-gma*(z-h)*eval(subs(k213*cos(delm213)));
else (y1>=(1-(1/(9*beta))) && y1<=1)
p2(i)=-gma*(z-h)*eval(subs(k201*cos(delm201)));
end
end
else
p2(i)=-k3*cos(delm3)*gma*(z-h-h1)-gma*k3*h1*cos(delta);
for y2=0:0.01:1
if(y2>=0 && y2<=(1/(9*alfa)))
p1(i)=gma*(z-h-h1)*eval(subs(k401*cos(delm401)))+(gma*(h1+h)+q)*k4R0*cos(delta/2);
elseif(y2>=(1/(9*alfa)) && y2<=(1/alfa))
p1(i)=gma*(z-h-h1)*eval(subs(k413*cos(delm413)))+(gma*(h1+h)+q)*k4R0*cos(delta/2);
else(y2>=(1/alfa) && y2<=1)
p1(i)=gma*(z-h-h1)*k43*cos(delm43)+(gma*(h1+h)+q)*k4R0*cos(delta/2);
end
end
end
end
z=0:0.1:L;
% subplot(2,1,1);
plot(p1,z,p2,z)
grid on
set(gca, 'YDir','reverse')

Raj on 4 Dec 2019
"The code is taking only last values, i.e., y1 and y2=1" - No it is not taking only the last values. The for loop for both y1 and y2 runs for all values from 0 to 1 with increment of 0.1
Lets see one loop:
for y1=0:0.01:1
if (y1>=0 && y1<=(1-(1/beta)))
p2(i)=-gma*(z-h)*k23*cos(delm23);
elseif (y1>=(1-(1/beta)) && y1<=(1-(1/(9*beta))))
p2(i)=-gma*(z-h)*eval(subs(k213*cos(delm213)));
else (y1>=(1-(1/(9*beta))) && y1<=1)
p2(i)=-gma*(z-h)*eval(subs(k201*cos(delm201)));
end
end
The problem here is that you are overwriting the P2 (i) value for each value of y1 and what gets finally saved is the last computed value with y1=1. If you want to save all the P2 values for every value of y1 you will have to index it properly by using a 2D matrix.

Show 3 older comments
Raj on 5 Dec 2019
If you have updated the code as above, p1 and p2 are 2D matrices where as Z is a 1D array. You cannot use
plot(p1,z,p2,z)
directly. Obviously the dimensions will not match. You have to modify the code to extract the data from p1 and p2 which you want to plot and then use plot command.
Akshay Pratap Singh on 5 Dec 2019
Thanks @Raj. however may you please help me by writing the code.