## 135 Reputation

13 years, 232 days

so when we write

round(rand())

Thanks Carl

## matlab code...

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);

## maple code...

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

## MC method...

in terms of MC method

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

## yes it is...

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

## minimum...

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

## Thanks...

Many many thanks..

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

Regards

## plz check...

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?

## maple...

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

## @Preben Alsholm  Yes I just realize...

Yes I just realize it now it is complex

Many thanks

## @Preben Alsholm  Thanks for ur resp...

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

## w(t) not oscillatory?...

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??

 5 6 7 8 9 10 11 Page 7 of 12
﻿