mix symbolic with function handle

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;

Answer 1

[Walter Roberson]

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]

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]

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]

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.