MATLAB Answers

mix symbolic with function handle

9 views (last 30 days)
Simo
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;

  0 Comments

Sign in to comment.

Accepted Answer

Walter Roberson
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

  3 Comments

Simo
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
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
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.

Sign in to comment.

More Answers (0)


Translated by