piecewise for use in anonymous functions

17 views (last 30 days)
Pooyan Jafari on 5 Aug 2020
Commented: Walter Roberson on 9 Aug 2020
Hello,
clc,clear all,close all;
opengl('save', 'software')
syms beta theta r ;
B = 10;
H = 10;
C = 10;
L1 = 4;
Ri = B/2 + L1;
Rf = C + B/2 + L1;
Rbeta0 = Ri;
theta1 = atan((H/2)/(B/2+L1));
theta2 = pi - atan((H/2)/(B/2-L1));
Rt0 = piecewise((0<=theta) & (theta<theta1),(B/2+L1)/cos(theta),(theta1<=theta) & (theta<theta2),(H/2)/sin(theta),(theta2<=theta) & (theta<=pi),(B/2-L1)/cos(pi-theta));
Rmaxtb0 = (Rt0*Rbeta0)/Ri;
Rmaxtb00 = matlabFunction(Rmaxtb0)
Vmbeta0 = (2*(pi^2)*((B+(2*L1))^2)) / ((((2*L1*pi)+(pi*B))^2));
Vbeta0 = Vmbeta0 * (1 - ((r^2) / (Rmaxtb0^2)));
Vbeta00 = Vbeta0 * r;
Vbeta000 = matlabFunction(Vbeta00,'Vars',[theta r]);
Wsigmat = 2 * integral2(Vbeta000,0,pi,0,Rmaxtb00,'AbsTol',1e-5, 'RelTol',1e-5)
Error using symengine
Unable to generate code for piecewise for use in anonymous functions.
As you can see in my piecewise function there is 'theta' in the limits which is defined as 'syms' and is not defined prior
and I think it couldn't be written as if/else!! (I don't know if it can)
Is there any way to define a piecewise function as an anonymous function?
or any way to use an alternative for piecewise function??
thanks in advance for any help

Walter Roberson on 6 Aug 2020
Rmaxtb00 = matlabFunction(Rmaxtb0, 'File', 'Rmaxtb00.m', 'optimize', false);
The generated function in the .m file will use if/elseif to code the piecewise()
However, the generated if/elseif code will not be vectorized, so when you use it with integral() you need to use 'ArrayValued', true so that integral() does not pass in a vector of values to be evaluated.

Pooyan Jafari on 9 Aug 2020
clc,clear,close all;
syms beta theta r ;
D = 10;
C = 10;
L1 = 4;
Ri = D/2 + L1;
Rf = C + D/2 + L1;
Rbeta = Ri + (2*beta*(Rf-Ri)/pi);
Rt0 = (2*L1*cos(theta) + sqrt((D^2) - (2*(L1^2)) + (2*(L1^2)*cos(2*theta))) )/2;
Rmaxtb = (Rt0*Rbeta)/Ri;
Rmaxtbb = matlabFunction(Rmaxtb);
Vmbeta = (2*(pi^2)*((D+(2*L1))^2)) / ((((2*L1*pi)+(pi*D)+(4*beta*C))^2));
Vmbeta0 = (2*(pi^2)*((D+(2*L1))^2)) / ((((2*L1*pi)+(pi*D))^2));
Vmbetapi2 = (2*(pi^2)*((D+(2*L1))^2)) / ((((2*L1*pi)+(pi*D)+(4*(pi/2)*C))^2));
Vbeta = Vmbeta * (1 - ((r^2) / (Rmaxtb^2)));
a_coeff = 512*r*C*((D+(2*L1))^2);
b_coeff = (L1*cos(theta)) * ((beta*C + 0.25*pi*D + 0.5*L1*pi)^2) * sqrt((D^2) - (2*(L1^2)) + (2*(L1^2)*cos(2*theta)));
c_coeff = (L1^2) * ((beta*C + 0.25*pi*D + 0.5*L1*pi)^2) * cos(2*theta);
d_coeff = ( ((pi*L1*(D-(2*r))) / 2) + (( (((0.25*D) - (0.5*r))*pi) + (beta*C)) * D) ) * ( ((pi*L1*(D+(2*r))) / 2) + (( (((0.25*D) + (0.5*r))*pi) + (beta*C)) * D) );
e_coeff = (((2*L1*pi)+(pi*D)+(4*beta*C))^5) * (r*cos(theta) - Rf);
f_coeff = (4*(L1^2)*cos(2*theta)) + (4*L1*cos(theta)*sqrt((D^2) - (2*(L1^2)) + (2*(L1^2)*cos(2*theta)))) + (D^2);
Vr = (-1 * a_coeff * (b_coeff+c_coeff+(d_coeff/4)) * (pi^2)) / (e_coeff*f_coeff);
epsrr = diff(Vr,r);
epsrtheta = (1/(2*r)) * (diff(Vr,theta));
epsrbeta = (1/(2*(Rf-(r*cos(theta))))) * ( (Vbeta*cos(theta)) + (diff(Vr,beta)) + ( (Rf-(r*cos(theta)))*(diff(Vbeta,r)) ) );
epsthetatheta = Vr/r;
epsthetabeta = (1/(2*r*(Rf-(r*cos(theta))))) * ( (-1*Vbeta*r*sin(theta)) + ( (Rf-(r*cos(theta)))*(diff(Vbeta,theta)) ) );
epsbetabeta = ( 1/(Rf-(r*cos(theta))) ) * ( (-1*Vr*cos(theta)) + (diff(Vbeta,beta)) );
g_coeff = (epsrtheta^2) + (epsthetabeta^2) + (epsrbeta^2);
h_coeff = (epsrr * epsthetatheta) + (epsthetatheta * epsbetabeta) + (epsbetabeta * epsrr);
p = g_coeff - h_coeff;
i_coeff = (epsrr * epsthetatheta * epsbetabeta) + (2*epsrtheta*epsthetabeta*epsrbeta);
j_coeff = (epsbetabeta * (epsrtheta^2)) + (epsrr * (epsthetabeta^2)) + (epsthetatheta * (epsrbeta^2));
q = i_coeff - j_coeff;
epsmax = 2 * sqrt(p/3) * cos((1/3) * acos( (3*sqrt(3)*abs(q)) / (2*sqrt(p^3)) ));
epsmaxx = (2*epsmax) *r*(Rf-(r*cos(theta)));
epsmaxxx = matlabFunction(epsmaxx,'Vars',[beta theta r]);
temp = epsmaxxx(sym('r'),0,0);
limit(temp, r, 0, 'left')
limit(temp, r, 0, 'right')
Dcu = integral3(epsmaxxx,0,pi/2,0,2*pi,0,Rmaxtbb,'AbsTol',1e-5, 'RelTol',1e-5)
the result:
ans =
NaN
ans =
NaN
Dcu =
768.8840
Walter Roberson on 9 Aug 2020
test = diff(Vbeta,theta)
test =
piecewise(theta in Dom::Interval(0, 6257646083598987/9007199254740992), etc
Notice that the Interval given starts with 0. We are looking at MuPad internal notation here, and in that notation there is a difference between the scalar 0 and the "list" [0] . MATLAB treats 0 and [0] the same (except that the [] involves an extra call). However MuPad does not treat them the same: 0 is DOM::Numeric and [0] is DOM::List which is a data structure. In the context of DOM:Interval, an endpoint that numeric is Open Interval, and an endpoint that is DOM::List is closed on that side.
Thus the closed interval endpoint in Vbeta includes 0 and you can subs(Vbeta, theta, 0) and get something useful. However when you differentiate the fact that Vbeta is not defined before 0 means that the derivative at 0 exactly is undefined. When you
subs(diff(Vbeta, theta), theta, 0)
then you get nan because 0 is outside any interval set up by the piecewise that remains after the derivative.
since the denominator cos(theta) is nonzero at 0, then it is not required that we introduce a singularity at 0 in the derivative; we just extend the piecewise slightly so that the derivative exists in theory at 0, and proceed knowing that we will be able to evaluate the derivative at 0 exactly in the piecewise without triggering nan.
Walter Roberson on 9 Aug 2020
"I think because Rmaxtb is in the denominator of Vbeta, there shouldn't be any problem regarding differentiation (I am not sure);"
You are missing that when you calculate a derivative then for any given answer to considered correct, it must be valid on every compact infinitesimal interval near the point. However when you have a function with a closed boundary and define the value outside the interval to be nan (because that is what piecewise does), then although the derivative might be correct on "one side" of the boundary, on the other side of the boundary the derivative is undefined. There is no non-zero infinitesimal that you can go negative to, and at which a hypothetical derivative is defined.
Can we potentially match endpoints between two piecewise that are defined on different sides of the boundary, piecewise one way and piecewise the other way agreeing on a value where they meet? That is a valid approach for some cases, but this case says that to the left (negative) the function value is undefined (nan). No possible meeting of functions.
subs theta=0 into the diff(Vbeta, theta) and you will get nan, for the reasons described. So you cannot integrate the triple integral on any interval that includes theta = 0.
... not unless you extend the piecewise in a way that is consistent, such as just moving the lower boundary from zero to something negative but greater than π/2