How to convert symbolic piecewise expression into a function handle

I require to convert a symbolic piecewise expression to a function handle in matlab.

I tried using matlabFunction() but it throws a symengine error of unbalanced or unexpected parenthesis or bracket.

Any help will be highly appreciated.

[Walter Roberson]

Please show your attempt to use matlabFunction .

[Vishal Tripathi]

clc; syms T; syms E; syms v; syms z;
assume (v>0 & v<6000);
fun = int(exp(-E/(8.314*T)),T,300,T);
s = 18 * 10^3; %Mean
fun1 = (1/(v*sqrt(2*pi)))* exp(- ((E-s)^2)/(2*v^2));
a = 1.95e+13; %Frequency factor
b = 10/60; %Heating Rate
fun2 = int(exp((-a/b)*fun)*fun1,E,s-3*v,s+3*v);
alpha1 =(z - fun2);%Calculated Alpha
T=[350 400];
alphaexp = [0.01134 0.04317];
alpha2 =(subs(alpha1))
minfunc = (alpha2-alphaexp).^2
error1 = (sum(minfunc))
error2 = matlabFunction(error1)
error = @(x)error2(x(1),x(2))
[x,fval]=fminsearch(error,[2000 1])

While calling the matlabFunction, the error shown is

Error using symengine>makeFhandle (line 109) Error: Unbalanced or unexpected parenthesis or bracket.

Error in symengine (line 60)

Error in sym/matlabFunction (line 125) g = symengine('makeFhandle',varnames,body);

Error in i2 (line 15) error = matlabFunction(error1)

Sorry for the late reply, and any help is badly needed.

Thanks.

P.S. - EDIT - The code is edited as per the discussions below.

[Walter Roberson]

In

alpha2 =(subs(alpha1))

which variable are you trying to substitute, and what value are you trying to substitute for it? You have not created any numeric variables with the same name as symbolic variables so there would not be anything obvious for the subs() to do. I recommend always be specific about what you are substituting for and which new value is being substituted in.

I do not see any piecewise definition in what you posted.

[Vishal Tripathi]

Since I have defined the values for 'T' just before the subs statement, I expect it to be substituted in alpha1.

After the integration, the function alpha1 has T, v and z as variables. Out of which I am substituting the values of T. Thus wanting to minimize the 'minfunc' with the minimum values of v and z.

Coming to the piecewise bit, the integration yields a piecewise function. Thus the problem.

P.S. - Since I do not receive notifications for the comments, I tend to reply late. Apologies for the inconvenience.

[Walter Roberson]

There is no piecewise if you add the assumption that T >= 300

I recommend you use a different name for the upper bound, such as

syms Tmax
assume(Tmax>300)
fun = int(exp(-E/(8.314*T)),T,300,Tmax);

Now by the time you get to error1, you have a function in which both z and v are free variables. Is it correct that you are trying to minimize over both variables simultaneously?

[Vishal Tripathi]

I think the integration yields a piecewise function because of E or v, and not because of T. However I think I may have found a solution to my problem. [by specifying a bracket for v, so that no piecewise function is formed]

And yes, I'm trying to minimize over both variables simultaneously.

Sadly another problem, which has popped in is that while using fminsearch the function considers E as an undefined function/variable, whereas I have clearly indicated in fun2 that the integration is over the variable E, with limits containing v. (so E should not even exist in the final error function)

Is there a problem with fminsearch, or am I doing some mistake somewhere ?

P.S - I have edited the code above for your reference.

Answer 1

[Walter Roberson]

If you write it out with the final integration step not done, and you expand() the terms and combine exponents and simplify, you can show that the minimum must occur when z is the mean of the alphaexp.

Near v = 0 and E = 18000, the integrals are very steep. The key point is near v = 0.002486920643673647643722704 where just a hair less gets you an integral that is so large that it effectively cannot be computed, but a hair more gets you an integral that is amazingly small -- for example values on the order of 1/10^(10^13). I think it means that the minimum is going to be at your upper bound, v < 6000 but the difference between that and any v > 0.002486920643673647643722704 would need to be computed to billions of digits to determine the difference. The inherent inaccuracy in your constants such as 8.314 do not justify those kinds of computations.

[Vishal Tripathi]

Apologies for the values. These were only decoys to go through the errors in the main code (which was a bit long)

I would however appreciate if you could explain the expand, combining and simplifying step. I am attaching the main code, with real values.

clc; clearvars; clear all;
syms T; syms E; syms v1; syms v2; syms v3;
assume(v1>0 & v1<50000);
assume(v2>0 & v2<50000);
assume(v3>0 & v3<60000);
b = 10/60;
fun = int(exp(-E/(8.314*T)),T,300,T);
s1 = 175.6 * 10^3;
fun11 = (1/(v1*sqrt(2*pi)))* exp(- ((E-s1)^2)/(2*v1^2));
a1 = 10^14.52;
fun12 = int(exp((-a1/b)*fun)*fun11,E,s1-3*v1,s1+3*v1);
alpha1 = 1 - fun12;
s2 = 185.4 * 10^3;
fun21 = (1/(v2*sqrt(2*pi)))* exp(- ((E-s2)^2)/(2*v2^2));
a2 = 10^13.64;
fun22 = int(exp((-a2/b)*fun)*fun21,E,s2-3*v2,s2+3*v2);
alpha2 = 1 - fun22;
s3 = 195.4 * 10^3;
fun31 = (1/(v3*sqrt(2*pi)))* exp(- ((E-s3)^2)/(2*v3^2));
a3 = 10^13.98;
fun32 = int(exp((-a3/b)*fun)*fun31,E,s3-3*v3,s3+3*v3);
alpha3 = 1 - fun32;
alpha = (alpha1 + alpha2 + alpha3)/3
alphaexp=[0.01134 0.04317];% 0.06494 0.08783 0.17053 0.32533 0.49142 0.55575 0.59242 0.6367 0.678 0.71621 0.75124 0.78442 0.81727];
T = [350 400]; %T = [350:50:1050];
minfunc = (subs(alpha)-alphaexp).^2
error1 = sum(minfunc)
error = matlabFunction(error1)
[xfinal,fval] = fminsearch(@(x)error(x(1),x(2),x(3)),[4300 3500 32000])
Variance = xfinal
Error = sqrt(fval)

[Walter Roberson]

Let me put this more firmly: your system is not doing the integration because your formula is too complex to integrate. It also cannot do effective numeric integration because your slope is too steep. You need to go back and reason about what is happening, like I did, rather than trying for an automatic solution that is not practical to achieve.

[Vishal Tripathi]

I understand what you are trying to say. However, what I cannot understand is that if I put the values of v1, v2 and v3, the integration is happening, and I am getting smooth results. Same is the case if I use for numerical integration (using vpa and defining v1, v2 and v3 specifically). Instead when I use fminsearch to optimise for same values, the variable 'E' is being treated as an alien function/variable which it should not according to me. This is what I am unable to grasp and need help for.

I am very sorry for being so ignorant, but since I am very new to MATLAB I request a bit of patience.

P.S - Kindly refer to the above code for computation.

[Walter Roberson]

Please do not name a variable "error" as that interferes with using the critical MATLAB function error() and makes it difficult for other people to read the code.

Please show the output of error1 and also the output of the matlabFunction call.

I have a vague memory that there might be a bug in the conversion of int() to function handle.