@Carl Love 

Thanks for ur fast reply..

so when we write

round(rand())

@Carl Love 

Thanks Carl

@YasH 

this is the matlab is it possible to rewrite it in simple maple code

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

B = 0;

M = []; % The total magnetic field of the system

E = []; % The total energy of the system

randTol = 0.1; % The tolerance, dampens the spin flip process

% First we generate a random initial configuration

spin = (-1).^(round(rand(N)));

% Then we let the system evolve for a fixed number of steps

for i=1:1000,

% Calculating the total spin of neighbouring cells

neighbours = circshift(spin, [ 0 1]) + ...

circshift(spin, [ 0 -1]) + ...

circshift(spin, [ 1 0]) + ...

circshift(spin, [-1 0]);

% Calculate the change in energy of flipping a spin

DeltaE = 2 * (J*(spin .* neighbours) + B*spin);

% Calculate the transition probabilities

p_trans = exp(-DeltaE);

% Decide which transitions will occur

transitions = (rand(N) < p_trans ).*(rand(N) < randTol) * -2 + 1;

% Perform the transitions

spin = spin .* transitions;

% Sum up our variables of interest

M = sum(sum(spin));

E = -sum(sum(DeltaE))/2; % Divide by two because of double counting

% Display the current state of the system (optional)

image((spin+1)*128);

xlabel(sprintf('J = %0.2f, M = %0.2f, E = %0.2f', J, M/N^2, E/N^2));

set(gca,'YTickLabel',[],'XTickLabel',[]);

axis square; colormap bone; drawnow;

end

% Count the number of clusters of 'spin up' states

[L, num] = bwlabel(spin == 1, 4);

I have the same problem have u rewrite it in maple??

in terms of MC method

First: choose arbitrary configration of spoin (up or down)

how can I start with this line

@Carl Love 

is there in maple a build in packege of monte calro in statistical mechanics?

@Markiyan Hirnyk @Kitonum

Actually the expression it is not sin(), 

plot3d(Re(f),Delta=-10..10,lambda=0..1, lightmodel=light2,axes=boxed,title=tit,view=0..0.005);

minimize(Re(f), Delta=-10..10, lambda=0..1, location);
Warning, computation interrupted

is there any way to get the minimum value from the plot or any comand can help other than minimze(this took time without any results)

the expression of Re(f) is too lengthy to write it here.

Regards

@Carl Love 

Many many thanks.. 

This is the method of monecalo.. are you familar how to apply it to any function.

Regards

Carl Love 9881

please wouldu plz advise me

@gkokovidis 

thanks for ur reply

size:=40;
40
limit:=size+1;
Error, attempting to assign to `limit` which is protected
x:=1;
ud:=0;
dd:=0;

while(x<(limit))
line([0, size], [x, x]);
x:=x+1;
end;

Is this right?

@gkokovidis 

How can I made the changes of the above code to be aplicable in Maple

@Preben Alsholm 

Yes I just realize it now it is complex

Many thanks

@Preben Alsholm 

Thanks for ur respond.

I think the above lim(sz,t goes to infinity) it isnot accurate in mathematics so I try to do it numericaly as this:

Delta:=5:omega:=10^6:alpha:=M/(2*omega):
N:=1:N1:=1+2*N;M:=sqrt(N*(N+1));

dsys :={diff(x(t),t)=-(N1-I*Delta-2*M*exp(-2*I*omega*t))*x(t), diff(y(t),t)=-(N1+I*Delta-2*M*exp(2*I*omega*t))*y(t), diff(z(t),t)=-(N1+M*cos(2*omega*t))*z(t)-1}:

eqs:=subs(x(t)=x,y(t)=y,z(t)=z,rhs~(dsys))=~0;
res:=solve(eqs,{x,y,z});
eqs2:=subs(x(t)=x1+I*x2,y(t)=y1+I*y2,z(t)=z,rhs~(dsys))=~0;
eqs2R:=(evalc@Re)~(eqs2);
eqs2I:=(evalc@Im)~(eqs2);
solve(eqs2R union eqs2I,{x1,x2,y1,y2,z});
eval(eqs2,%);
simplify(%) assuming real;

subs(RootOf(conjugate(_Z)+_Z)=0+I*x2,res);

Here I got the exact value of z(t).

Now my question is how to use z(t) to calculate: z(t)*w1 and plot it.

where:

I3:=Sum((BesselJ(n,-2*alpha)*(-1)^n+BesselJ(n,-2*I*alpha))*exp(-(N1-I*Delta+2*n*omega)*t),n=-infinity..infinity):
w:=inttrans[laplace](I3,t,s):
w1:=subs(s=I*d,w):

Many Many thanks

@Preben Alsholm 

Many many thanks..

When I run this code only w(t) not oscillatory ... which is not true because of the precence of the terms cos() and sin()

Digits:=15:
x(0):=-1:y(0):=0:z(0):=conjugate(y(0)):N:=1:Delta:=5: N1:=1+2*N:M:=sqrt(N*(N+1)):
ini1:=u(0)=Re(y(0)), v(0)=Im(z(0)),w(0)=x(0);
omega:=10^(6):
dsys1 :=diff(w(t),t)=-(N1+M*cos(2*omega*t))*w(t)-1+2*u(t)*cos(2*omega*t)+2*v(t)*sin(2*omega*t), diff(u(t),t)=-N1*u(t)+Delta*v(t)+(2*M*u(t)-N1-w(t))*cos(2*omega*t)-2*M*v(t)*sin(2*omega*t), diff(v(t),t)=-N1*v(t)-Delta*u(t)+(2*M*u(t)-N1-w(t))*sin(2*omega*t)+2*M*v(t)*cos(2*omega*t);

dsol1 :=dsolve({dsys1,ini1},numeric,abserr=1e-9, relerr=1e-8,maxfun=0);
plots:-odeplot(dsol1,[[t,w(t)]],0.9..1,axes=boxed,tickmarks=[6, 6],axes=boxed,titlefont=[SYMBOL,12]);

If u run the above would be the same??