# mix symbolic with function handle

9 views (last 30 days)
Simo on 24 Jun 2020
Commented: Walter Roberson on 24 Jun 2020
hi
I've to write a function with parameter b, where an integral is in this function. how can i do it? thanks for any advice
L=2; N=L*L; j=1;
Kb=1;
syms b
k=2*sinh( 2 * b * j )/(cosh( 2 * b * j))^2;
intg=@(x)(log ( 0.5 * ( 1 + sqrt(1 - (k(b)).^2 .* (sin(x))^2 ))));
INTEGRAL=(1/(2*pi)) * integral(intg,0,pi);
Z(b)=(2* cosh(2 * b * j) * exp(INTEGRAL))^N;
E_form(b)=-diff(log(Z))/N;

Walter Roberson on 24 Jun 2020
integral() can never be used with an unresolved symbolic variable.
If you have unresolved symbolic variables then you need to use int()
Also note that you created k as an expression rather than a function, so you should be using just k instead of k(b)
L = sym(2);
N = L*L;
j = sym(1);
Kb = sym(1);
syms b x
k = sym(2) * sinh( sym(2) * b * j )/(cosh( sym(2) * b * j))^sym(2);
intg = @(x)(log ( sym(0.5) * ( 1 + sqrt(1 - (k).^2 .* (sin(x))^2 ))));
INTEGRAL = (sym(1)/(sym(2)*sym(pi))) * int(intg, x, sym(0), sym(pi));
Z(b) = (sym(2)* cosh(sym(2) * b * j) * exp(INTEGRAL))^N;
E_form(b) = -diff(log(Z))/N;
improved_E_form = simplify(E_form);
You will notice that the result still has int() in it. Some (many) symbolic patterns are difficult find find closed form integrals for.
You were probably trying to use numeric integrals to avoid ending up with symbolic integrals, but that simply will not work when the expression has an unresolved variable in it: int() will solve whatever parts it can but what is left cannot be handled with numeric integration .
There is a closed form formula, by the way, one involving the Elliptic E and Elliptic K integrals:
if b == 0
result = 0;
else
Pi = pi;
result = -512*sinh(2*b)*(1/4*(cosh(b)^2-1/2*cosh(2*b)-1/2)^2*(cosh(b)^2+1/2*cosh(2*b)-1/2)^2 ...
* (cosh(2*b)^2-2)*((cosh(b)^4-cosh(b)^2)*cosh(2*b)^2-(cosh(b)^2-1/2)^2)^2*(cosh(b)^2-1/2)^4 ...
* EllipticE(4*cosh(b)/(2*cosh(b)^2-1)^2*sinh(b))+(cosh(b)+1)*cosh(b)^2 ...
* ((cosh(2*b)^2-2)*(cosh(2*b)-1)*(cosh(2*b)+1)*(cosh(b)^2-1/2)^8 ...
* EllipticK(4*cosh(b)/(2*cosh(b)^2-1)^2*sinh(b))+cosh(2*b)^4*(cosh(b)+1) ...
* (-1/128*cosh(2*b)^6+1/64*cosh(2*b)^4+(cosh(b)+1)*cosh(b)^2*(cosh(b)-1)*(cosh(b)^2-1/2)^4) ...
* cosh(b)^2*(cosh(b)-1)*Pi)*(cosh(b)^4-cosh(b)^2-1/4)^2*(cosh(b)-1)) / Pi ...
/ sinh(b)^6 / cosh(b)^6 / (4*cosh(b)^4-4*cosh(b)^2-1)^2/(2*cosh(b)^2-1)^9;
end

Simo on 24 Jun 2020
it works, thanks
I don't understand: is there a way to speed the code?
Cv_form(b)=-diff(E_form) * Kb * b^2/N;
for iT=1:50
T(iT)=0.1*iT;Kb=1; KbT=Kb*T(end); beta=1/KbT;
E_th(iT)=E_form(beta);
Cv_th(iT)=Cv_form(beta)/(Kb*T(iT)*T(iT));
end
Walter Roberson on 24 Jun 2020
You are expanding T dynamically. You should be using vectorized calculations.
iT = 1 : 50;
T = 0.1 * iT;
Kb = 1;
KbT = Kb * T;
beta = 1 ./ KbT;
E_th = E_form(beta);
Cv_th = Cv_form(beta) ./ (Kb .* T.^2);
Walter Roberson on 24 Jun 2020
You could also experiment with using vpaintegral() instead of int() . I am not certain at the moment how that will interact with your diff() calls.