Definite Integrals in for-loop

I'm trying to setup up this loop. But before I do that I need to declare the functions:

n_grid   =    20;  %  Size  of    rid
Ms       =    [];
Js       =    [];
Ns       =    [];
Es       =    [];
J        =    linspace(0,1,1000);
kappa = 2.*sinh(2*J)./(cosh(2*J)).^2;
kappam = 2*(tanh(2*J)).^2-1;
z = exp(-2*J);
K1 = int(sqrt(1- kappam.^2 * sin(phi).^2), phi= 0.. pi/2)
% Calculating Energy & Magnetisation as a function of interaction strength
E = -J.*(coth(2*J)).*(1+(2/pi)*kappam.*K1(kappam));

An error for K1 shows as follows:

Undefined function 'K1' for input arguments of type 'char'.

Error in ising (line 12)
K1 := int(sqrt(1- kappam.^2 * sin(phi).^2), phi= 0.. pi/2)

What's wrong with my syntax?

Answer 1

[Andrei Bobrov]

Please try code below

J = linspace(0,1,1000);
kappam = 2*(tanh(2*J)).^2-1;

f = @(phi)sqrt(1- kappam.^2 * sin(phi).^2);
K1 = integral(f,0,pi/2,'ArrayValued',true);

[T]

Hi Andrei, thanks for your help. I receive the following error:

Subscript indices must either be real positive integers or logicals.

Subscript indices must either be real positive integers or logicals.
Error in ising (line 14)
E = -J.*(coth(2*J)).*(1+(2/pi)*kappam.*K1(kappam));

[Andrei Bobrov]

What is it K1(kappam), may by just K1?

[T]

This code is from page 11 of this document: https://notendur.hi.is/jeg1/Ising.pdf

K1 is an integral from page 3 (2.14)

[Andrei Bobrov]

Hi T!

two variants:

first:

f = @(phi)1./sqrt(1- kappam.^2 * sin(phi).^2);
K1 = integral(f,0,pi/2,'ArrayValued',true);
E = -J.*(coth(2*J)).*(1+(2/pi)*kappam.*K1);

second:

f = @(k,phi)1./sqrt(1- k.^2 * sin(phi).^2);
K1 = @(k)integral(@(phi)f(k,phi),0,pi/2,'ArrayValued',true);
E = -J.*(coth(2*J)).*(1+(2/pi)*kappam.*K1(kappam));

[T]

In using the second part:

Error using integralParseArgs (line 12)
Argument 'Ar' did not match any valid parameter of the parser.

Error in integral (line 88)
opstruct = integralParseArgs(varargin{:});

Error in @(k)integral(@(phi)f(k,phi),0,pi/2,'Ar',1)

Error in ising (line 14)
E          =          -J.*(coth(2*J)).*(1+(2/pi)*kappam.*K1(kappam));

Nevertheless, I think 'k' should actually be kappa as defined in 2.12 on page 3 and in the code:

J          =          linspace(0,1,1000);
kappa          =          2.*sinh(2*J)./(cosh(2*J)).^2;

[Andrei Bobrov]

corrected

first:

f = @(phi)1./sqrt(1- kappa.^2 * sin(phi).^2);
K1 = integral(f,0,pi/2,'ArrayValued',true);
E = -J.*(coth(2*J)).*(1+(2/pi)*kappam.*K1);

second:

f = @(k,phi)1./sqrt(1- k.^2 * sin(phi).^2);
K1 = @(k)integral(@(phi)f(k,phi),0,pi/2,'ArrayValued',true);
E = -J.*(coth(2*J)).*(1+(2/pi)*kappam.*K1(kappa));

[T]

That works:

n_grid   =          20;          %  Size  of    rid
Ms          =          [];
Js          =          [];
Ns          =          [];
Es          =          [];

J          =          linspace(0,1,1000);
kappa          =          2.*sinh(2*J)./(cosh(2*J)).^2;
kappam          =          2*(tanh(2*J)).^2-1;
z          =          exp(-2*J);

f = @(kappa,phi)1./sqrt(1- kappa.^2 * sin(phi).^2);
K1 = @(k)integral(@(phi)f(kappa,phi),0,pi/2,'ArrayValued',true);
%          Calculating          Energy          &          Magnetisation          as          a          function          of          interaction          strength
E          =          -J.*(coth(2*J)).*(1+(2/pi)*kappam.*K1(kappam));
M1          =          ((1+z.^2).^(1/4).*(1-6*z.^2+z.^4).^(1/8))./(1-z.^2).^(1/2);
M2          =          -((1+z.^2).^(1/4).*(1-6*z.^2+z.^4).^(1/8))./(1-z.^2).^(1/2);
M1(1,1:441)          =          0;          %When          J          <          J_c,          M          =          0
M1(1,1:441)          =          0;          %When          J          <          J_c,          M          =          0

for          i          =          1:12000

%while          (1),
%          Choose          a          random          value          between          0          and          1          for          the          interaction          strength
J          =          rand()+1e-10;
%          Perform          a          simulation
[M,          N,          E]          =          ising2(n_grid,          J);
%          Records          the          results
Ms          =          [Ms          M/(n_grid^2)];
Es          =          [Es          E/(n_grid^2)];
Ns          =          [Ns          N];
Js          =          [Js          J];
i          =          i+1;
end

Undefined function 'ising2' for input arguments of type 'double'.

Error in ising (line 27) [M, N, E] = ising2(n_grid, J);

[Walter Roberson]

If you look on page 12 of that PDF you will see that even they have an error message about a missing ising2 function or variable!

Your file that you have named "ising.m" should instead be named "montecarlo.m" . You should be taking the ising.m from that pdf and storing that as ising.m . And then you should create

function [M, num, E] = ising2(N, J)
  [M, num, E] = ising(N, J);

and store that as ising2.m

[T]

Done. Three errors arise:

Error: File: ising.m Line: 12 Column: 1
Function definitions are not permitted in this context.

Error in ising2 (line 2)
  [M, num, E] = ising(N, J);

Error in montecarlo (line 27)
[M,          N,          E]          =          ising2(n_grid,          J);

Were you able to reproduce the plots?

[Walter Roberson]

I had no problems with the software once I had done what I indicated above needed to be done.

You might notice that the code in ising.m is spaced oddly. That is because I copied it directly from the pdf and did not bother to change the spacing.

[T]

Strange. Those files work thanks.