piecewise for use in anonymous functions

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

Answer 1

[Walter Roberson]

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]

thanks for your answer

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, 'File', 'Rmaxtb00.m', 'optimize', false)

Vmbeta0 = (2*(pi^2)*((B+(2*L1))^2)) / ((((2*L1*pi)+(pi*B))^2));
Vbeta0 = Vmbeta0 * (1 - ((r^2) / (Rmaxtb0^2))) * r;
Vbeta00 = matlabFunction(Vbeta0, 'File', 'Vbeta00.m', 'optimize', false,'Vars',[theta r])

Wsigmat = 2 * integral2(Vbeta00,0,pi,0,Rmaxtb00,'ArrayValued',true,'AbsTol',1e-5, 'RelTol',1e-5)

I got this error:

Error using integral2ParseArgs (line 13)

'ArrayValued' is not a recognized parameter. For a list of valid name-value pair arguments, see the documentation for this function.

[Walter Roberson]

ArrayValued is for integral() not integral2()

[Pooyan Jafari]

Thank you so much

I have converted it as two nested integral and it is working for my double integral.

I have two triple integral in the rest of my code;

the first one is easier and I used the method you have proposed and it is working almost fine

the second triple integral in my code is more complex than previous one and I got this error message:

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.1e+02. The integral may
not exist, or it may be difficult to approximate numerically to the requested accuracy. 

vpaintegral() is working for my problem and gives the correct answer, but the computation time is very high and is not practical for me!!!

would you please take a look at the code and see if you can find a way to solve the problem...

(In an article someone has done this with Matlab software with a more complex piecewise function than me but I don't know how...)

I would be grateful if you could help me in this matter

opengl('save', 'software')
syms beta theta r ;
B = 10;
H = 10;
C = 10;
L1 = 1;

Ri = B/2 + L1;
Rf = C + B/2 + L1;
Rbeta = Ri + (2*beta*(Rf-Ri)/pi);
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));

Rmaxtb = (Rt0*Rbeta)/Ri;
Rmaxtbb = matlabFunction(Rmaxtb, 'File', 'Rmaxtbb.m', 'optimize', false,'Vars',[theta beta]);

Vmbeta = (2*(pi^2)*((B+(2*L1))^2)) / ((((2*L1*pi)+(pi*B)+(4*beta*C))^2));
Vbeta = Vmbeta * (1 - ((r^2) / (Rmaxtb^2)));

a_coeff = 1 / ( (Rmaxtb^3) * ((int(Rmaxtb^2,theta,0,pi))^2) * r * (r*cos(theta)-Rf));
b_coeff = (((r^2) * Rmaxtb) - (Rmaxtb^3)) * (int(Rmaxtb * diff(Rmaxtb,beta),theta,0,pi/2));
c_coeff = ( int(Rmaxtb^2,theta,0,pi) ) * ( diff(Rmaxtb,beta) ) * ((r^2) + (Rmaxtb^2));
d_coeff = (H*(r-Rmaxtb)) * (B*(r+Rmaxtb)) * (b_coeff + c_coeff);
Vr = a_coeff * d_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, 'File', 'epsmaxxx.m', 'optimize', false,'Vars',[r theta beta]);

int1 = @(theta,beta) integral(@(r)epsmaxxx(r,theta,beta),0,Rmaxtbb(theta,beta), 'ArrayValued',1);
int2 = @(beta) integral(@(theta)int1(theta,beta),0,pi, 'ArrayValued',1); 
Dc = 2 * integral(int2, 0,pi/2, 'ArrayValued',1,'AbsTol',1e-5, 'RelTol',1e-5)

[Walter Roberson]

The most immediate problem is that the way you have defined Rt0 with piecewise, when diff(Vbeta,theta) is taken, the derivative is not defined at 0, so you get a NaN being generated (because piecewise gives nan for any value out of range.)

That is something you can define away by changing to

Rt0 = piecewise((-pi/2 < 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));

This extends the first case across 0, allowing the derivative to exist at 0 exactly.

That done, you still have problems. The (extended) integral does not converge at r = 0, theta = 0, beta = 0.

temp = epsmaxxx(sym('r'),0,0);
limit(temp, r, 0, 'left')
limit(temp, r, 0, 'right')

MATLAB will not be able to find those two limits, but if you put the expression through another package such as Maple, you will find that the left limit is -infinity and the right limit is +infinity .

The message about reaching the maximum number of intervals is generated while integrating with respect to r, when beta and theta are both small and positive, such as

       theta: 1.71978371262361e-06
        beta: 8.59891856311806e-07

By that point, the function being integrated has reached about 2E+17 and getting larger rapidly as r gets smaller.

Because it is increasing in value so quickly, and the values are all positive at the time the problem is detected, we know that this is not just a removable discontinuity being dealt with, and that exact way that we extended Rt0 to be continuous across 0 is not the cause of the problem (because theta is never negative in the calculation.)

So, I don't think you can caclulate the integral -- not unless you put in non-zero boundary conditions.

[Pooyan Jafari]

Thank you for your answer

I think because Rmaxtb is in the denominator of Vbeta, there shouldn't be any problem regarding differentiation (I am not sure);

test = diff(Vbeta,theta)

test =
 
piecewise(theta in Dom::Interval(0, 6257646083598987/9007199254740992), (3125302502557517*r^2*cos(theta)*sin(theta))/(549755813888*(40*beta + 12*pi)^2*((20*beta)/pi + 6)^2), theta in Dom::Interval(6257646083598987/9007199254740992, 316031275249939/140737488355328), -(28127722523017653*r^2*cos(theta)*sin(theta))/(3435973836800*(40*beta + 12*pi)^2*((20*beta)/pi + 6)^2), theta in Dom::Interval(316031275249939/140737488355328, pi), (28127722523017653*r^2*cos(theta)*sin(theta))/(2199023255552*(40*beta + 12*pi)^2*((20*beta)/pi + 6)^2))
 

I haven't seen any NaN in the result.

As you can see in the following code, similarly MATLAB is not able to find the limits you have mentioned.

but it is calculating the integral and is giving the correct result.

here and here you guided me how to use integral() instead of vpaintegral()

The only difference of these two codes is in Rt0 and the rest is the same.

Rt0 in circular shape is defined with only one equation and everything is working fine.

but in square shape (my code), it is a piecewise function!!

As I mentioned before, vpaintegral() is giving the correct answer and is working fine (just high calculation time)

so I think my formulation doesn't have any problem

[Pooyan Jafari]

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]

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]

"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