## How to make a two dimensional ising model?...

I would like to make a Monte Carlo 2-D ising model but I am clueless on how I should start and so I would appreciate any advise/guidlines on this.

## Solution of definite integration...

```restart;
Digits := 15:
L := 1:
E := 100:
nu := 0.2:
G := E/2.6:
h := 0.1:
b := 0.1:

s := -E*(diff(w(x), x, x))*sinh(sqrt(2*Pi^2*N*(1+nu)/L^2)*y)/(sqrt(2*Pi^2*N*(1+nu)/L^2)*cosh(sqrt(Pi^2*N*(1+nu)/(2*L^2))*h)):

t := G*(diff(w(x), x))*(1-cosh(sqrt(2*Pi^2*N*(1+nu)/L^2)*y)/cosh(sqrt(Pi^2*N*(1+nu)/(2*L^2))*h)):

integrand := b*(int(t^2/(2*G)+s^2/(2*E), y = -(1/2)*h .. (1/2)*h))-(1/2)*E*b*h^3*evalf(Pi^2)*N*(diff(w(x), x))^2/(12*L^2):

integrand := subs(diff(w(x), x, x) = S, diff(w(x), x) = F, w(x) = Z, integrand):

EQ := subs(S = diff(w(x), x, x), F = diff(w(x), x), Z = w(x), diff(integrand, Z))-(diff(subs(S = diff(w(x), x, x), F = diff(w(x), x), Z = w(x), diff(integrand, F)), x))+diff(subs(S = diff(w(x), x, x), F = diff(w(x), x), Z = w(x), diff(integrand, S)), x, x):

W := rhs(dsolve(EQ)):

u1 := (int((E*sinh(sqrt(2*Pi^2*N*(1+nu)/L^2)*y)/(sqrt(2*Pi^2*N*(1+nu)/L^2)*cosh(sqrt(Pi^2*N*(1+nu)/(2*L^2))*h)))^2, y = -(1/2)*h .. (1/2)*h))*(int((diff(W, x, x))^2, x = 0 .. L))/(2*E):

u2 := (int(G^2*(1-cosh(sqrt(2*Pi^2*N*(1+nu)/L^2)*y)/cosh(sqrt(Pi^2*N*(1+nu)/(2*L^2))*h))^2, y = -(1/2)*h .. (1/2)*h))*(int((diff(W, x))^2, x = 0 .. L))/(2*G):

U := simplify(u1+u2-(1/2)*E*b*h^3*evalf(Pi^2)*N*(int((diff(W, x))^2, x = 0 .. L))/(12*L^2))```

It seems that solving of the above integration is a very time consuming process.

Please propose a way to solve above integration, if it is possible.

Thanks

## How do I ignore terms containing certain combinati...

I need to simplify a differential equation by ignoring higher order terms. The terms are of the following form:

180.495*theta_dot*beta_dot*zeta_dot*psi_dot*sin(beta)*cos(delta__3)*sin(delta__3)^3*zeta_dot*cos(zeta)*zeta*beta*sin(psi)*p

Here beta, zeta and theta are hinge deflections for a helicopter blade and thus can be assumed to be small. Everything other than delta__3 is a function of time, _dot represents a derivative wrt time. Higher order terms, such as those containing betan, zetan, thetan where n > 1, and terms containing combinations of the hinge deflections, i.e. betaj.zetak, betaj.thetak or zetaj.thetak where j,k >= 1 need to be neglected. Since beta, zeta and theta are functions of time, and the equation contains derivatives, maple functions like 'degree' and 'match' don't seem to work. Also, the solution suggested over here does not seem to be suitable for this problem, since my equation is very long (~70 pages) and the various derivatives appearing in it are not known.

## How to use an array/matrix as bcs for a pde equati...

First: Is it possible to solve an pde equation/system where bcs are an array or matrix?

Second: Is it possible to get a discrete response in an array/matrix for this pde equation/system?

As an example I applied the heat conduction equation to a bar of 50 cm in length. Initially the bar has a uniform temperature field of 20ºC. The ends are maintained at a temperature of 0°C over time. I fixed a unitary thermal diffusivity.

restart:
k:=1:
pde:=k*diff(v(x,t),x\$2)=diff(v(x,t),t):
bc:={v(x,0)=20,v(0,t)=0,v(50,t)=0}:
sol:=pdsolve(pde,bc,numeric,time=t,range=0..50):
p1:=sol:-plot(t=0,numpoints=50,color=red):
p2:=sol:-plot(t=20,numpoints=50,color=blue):
p3:=sol:-plot(t=50,numpoints=50,color=red):
p4:=sol:-plot(t=150,numpoints=50,color=blue):
p5:=sol:-plot(t=300,numpoints=50,color=red):
plots[display]({p1,p2,p3,p4,p5});

I will put the two questions in another way:
I would like to insert the boundary conditions not as algebraic functions, but rather as array/matrix. It is possible?
I would like to result in not a procedure, but rather as array/matrix. It is possible?

## Sum not working as intended....

I'm working with some sum values, but for some reason that I can't figure out, this sum always returns 1.
Have I made some mistake in the way it's supposed to be typed in or why does it return 1? Below the sum I've filled out the spots manually to show how I want them to be.

 >
 (1)
 >
 (2)
 >

## How to make 2nd derivative smoother?...

The issue concerns calculation of 2nd derivative of the numerical solution of 2nd order BVP for an ordinal DE.

Lets consider a test problem:

-y’’(t)+y(t)=t*sin(5*t)

y(0)=0.2, y(1)=1

Numerical solution was obtained by

`dsolve([-diff(F(x), x\$2)+F(x)= x*sin(5*x), F(0)=0.2, F(1)=1], [F(x)], type = numeric, 'output' = Array([seq(k/5, k=0..5)]));`

dsolve returns the values of the solution and its 1st derivative (for the test example, a solution can be obtained analytically, but for general case I require numerical solution)..

It is needed to calculate the 2nd derivative of the solution.

I tried to use the Bessel method, when the solution’s 2nd derivative is calculated as the 2nd derivative of 5th degree polynomial having in consequent points t1, t2, t3 the values y(t1), y(t2), y(t3) and 1st derivatives y’(t1), y’(t2), y’(t3), all obtained from the numerical solution. Later I use Hermit piecewise-polynomial interpolation, based on values of the solution and its 1st and 2nd derivatives (polynomial of 5th degree).

Unfortunately, 2nd derivative of such interpolation has a large non-smoothness.

Here, we see the solution. It is smooth.

Here, 1st derivative. Still smooth enough.

Here, 2nd derivative. Typical view with large and sharp teeth.

Maybe, there exists a simple method for calculation of more smooth 2nd derivative?

Also. it is desirable that it would be embedded in Maple.

## How I can solve these differential equations?...

How I can solve these time delay  differential equations?

Best

Doc191.pdf

ny.mw

 >
 (1)
 >
 (2)
 >
 (3)
 >
 (4)
 >
 (5)
 >
 (6)
 >

## Possible to disable the editable feature in Maple?...

Hi guys and girls,

The new Maple 2018 have a feature in the buttom called editable. If you remove that mark in the buttom of the screen, then maple doc can't be edited. I teach students who have english as a second language and some remove the mark in the editable box by mistake. Such that they can't do there hand-ins etc.

So is there anyway to disable the editable box in the buttom of the screen?

Fred

## How to plot Contour...

with(plots); R1 := .1; R0 := .1; m := .1; a := .1; Ha := .1; Nt := .1; Nb := .1; Pr := 6.2; Le := .6; Bi := 1; Ec := .1; k := 1; r := .1; A := 1; fcns := {C(y), T(y), U(y), W(y)}; sys := diff(U(y), `\$`(y, 2))+(R1*(diff(U(y), y))-2*R0*W(y))*exp(a*T(y))-a*(diff(U(y), y))*(diff(T(y), y))-Ha*(U(y)+m*W(y))*exp(a*T(y))/(m^2+1)-U(y)/k+A*exp(a*T(y)) = 0, diff(W(y), `\$`(y, 2))+(R1*(diff(W(y), y))+2*R0*U(y))*exp(a*T(y))-a*(diff(W(y), y))*(diff(T(y), y))-Ha*(W(y)-m*U(y))*exp(a*T(y))/(m^2+1)-W(y)/k = 0, diff(T(y), `\$`(y, 2))+R1*Pr*(diff(T(y), y))+Pr*Ec*exp(-a*T(y))*((diff(U(y), y))*(diff(U(y), y))+(diff(W(y), y))*(diff(W(y), y)))+Nt*(diff(T(y), y))*(diff(T(y), y))+Pr*Ec*(U(y)*U(y)+W(y)*W(y))*exp(-a*T(y))/k = 0, diff(C(y), `\$`(y, 2))+Pr*Le*R1*(diff(C(y), y))+Nt*(diff(C(y), `\$`(y, 2)))/Nb-r*C(y) = 0; bc := U(0) = 0, W(0) = 0, C(0) = 0, (D(T))(0) = Bi*(T(0)-1), U(1) = 0, W(1) = 0, C(1) = 1, T(1) = 0; L := [.5, 1.0, 1.5, 2.0]; AP := NULL; for k to 4 do R := dsolve(eval({bc, sys}, Ha = L[k]), fcns, type = numeric, method = bvp[midrich], AP); AP := approxsoln = R; p1u[k] := odeplot(R, [y, U(y)], 0 .. 1, numpoints = 100, labels = ["y", "U"], linestyle = dash, color = black) end do; display({p1u[1], p1u[2], p1u[3], p1u[4]})

## how to fix plot range to see the change in results...

for small value of k no change show in plot... the resulted value is so small and the system consider it zero and show no change in plot.. then how to get the change in plot to change in k... ??? plz suggest

 >
 >
 >
 (1)
 >
 >
 >
 >
 >
 >
 >
 >
 > sys := ([seq])(seq(eq[i, j], j = 0 .. N), i = 1 .. N):
 > nops(sys);
 > vars:=indets(sys) minus {k}:
 > nn := Matrix(N+1, N+1,(i, j)-> u[i-1, j-1]):
 > ##
 > p:=proc(kk) local u_res,A;
 > u_res:=solve(eval(sys,k=kk),vars);
 > A:=eval(nn,u_res);
 > plots:-matrixplot(A)
 > end proc;
 > ## Testing p for k=0.001:
 > p(0.001);
 > ## Animating the plot for k=0.0001..0.001:
 > plots:-animate(p,[k],k=0.0001..0.001);
 >
 >
 >

 >
 >
 >
 (1)
 >
 >
 >
 >
 >
 >
 >
 >
 > sys := ([seq])(seq(eq[i, j], j = 0 .. N), i = 1 .. N):
 > nops(sys);
 > vars:=indets(sys) minus {k}:
 > nn := Matrix(N+1, N+1,(i, j)-> u[i-1, j-1]):
 > ##
 > p:=proc(kk) local u_res,A;
 > u_res:=solve(eval(sys,k=kk),vars);
 > A:=eval(nn,u_res);
 > plots:-matrixplot(A)
 > end proc;
 > ## Testing p for k=0.001:
 > p(0.001);
 > ## Animating the plot for k=0.0001..0.001:
 > plots:-animate(p,[k],k=0.0001..0.001);
 >
 >
 >

## How do I solve a differential equation numerically...

Dear all, I had a challenge of using maple to obtain numerical solution of a volterra-integral equation using the method of reproducing kernel space, however I converted the integral equation to an ode and try to solve it numerically but obtatined an error message. Please how do I overcome it? the problem is as follows:

dsn := dsolve({diff(u(x), x) = e^x*ln(e)+(2/9)*e^3-1/9+int(y*u(y)^3, y = 0 .. 1), u(0) = 1}, numeric)

## Sum and with(SignalProcessing)...

Here in my code I am trying to differentiate a finite summation and the command Sum is giving me a simpler evaluation than the command sum.

But I have a problem, the command Sum doesn't work under signal processing tool. Is there away to make it work?

## Why does the code not work in Maple 2018.1?...

For some unknown reason, the code below does not work in Maple 2018.1, but works in Maple 2015 and Maple 2017 (the idea is taken from here

```restart;
with(plottools): with(plots):
V1,V2,V3,V4,V5,V6,V7,V8:=[0,-1,0],[0,0,0],[1,0,0],[1,-1,0],[0,-1,1],[0,0,1],[1,0,1],[1,-1,1]:  # The vertices of the cube
Faces:=[[V1,V4,V8,V5],[V5,V6,V7,V8],[V2,V3,V7,V6],[V1,V2,V3,V4],[V3,V4,V8,V7],[V1,V2,V6,V5]]: # The list of the faces
Colors:=[green, red,RGB(1, 0, 4),blue,grey,gold]: # The list of the colors
Cube[0]:=display([seq(polygon(Faces[i],color=Colors[i]),i=1..6)]):

for n from 1 to 7 do
F[n]:=t->rotate(Cube[n-1],t, [[0,n-1,0],[1,n-1,0]]):
Cube[n]:=rotate(Cube[n-1],-Pi/2, [[0,n-1,0],[1,n-1,0]]):
A[n]:=animate(display,[F[n](t)], t=0..-Pi/2,paraminfo=false);
od:

for m from 6 to 0 by -1 do
G[m]:=t->rotate(Cube[m+1],t, [[0,m,0],[1,m,0]]):
B[m]:=animate(display,[G[m](t)], t=0..Pi/2,paraminfo=false);
od:

C1:=display([seq(A[k], k=1..7)], insequence):
C2:=display([seq(B[k], k=6..0, -1)], insequence):
display([C1,C2], insequence, scaling=constrained, axes=normal);
```