## 70 Reputation

4 years, 234 days

## @Carl Love  Here's the correct...

Here's the corrected code:

restart:  with(Statistics):  Seed:= randomize(): # set up random numbers   N:= 1000: # number of samples   deltax:= 0.1: # max displacement nm   nm:= 10^(-9): # nanometres   kB:= 1.381*10^(-23)/nm^2: # Boltzmann const J nm-2 K-1   T:= 300.0: # temperature K   kBT:= kB*T: # initial thermal energy    V:= x-> 0.01*kBT*x^2 /2: # pot'l energy kg nm2 s&ndash;2   Etot:= 0.0: # initial energy <E>   E2tot:= 0.0: # initial <E2>  Xtot:=0.0:#initial <X>  X2tot:=0.0:#initial<X2>   x1:= 0.0: # first guess of x in nm   E1:= V(x1): # first guess of energy   for i from 1 to N do # start loop step (2)   x2(i):= x1 + (rand()/1e12 - 0.5)*deltax: # new x position   E2:= V(x2(i)): # new PE   DeltaE:= E2 - E1: # energy difference  # next line is Metropolis sampling   if DeltaE <= 0.0 or exp(-DeltaE/kBT) > (rand()/(1e12)) then   x1:= x2(i): # save new configuration   E1:= E2: # save new energy   end if:   Etot:= Etot + E1: # always add to total <E>  E2tot:= E2tot + E1^2: # add to total <E2>  Xtot:= Xtot + x1: # add to total <X>   X2tot:= X2tot + x1^2: # add to total <X2>   end do: # end loop step (3)   # convert result back to SI units  Rx:=AutoCorrelation(<seq(x2(i),i=1..N)>,1000): #compute autocorrelation of x(t)  RE:=AutoCorrelation(<seq(E2(i),i=1..N)>,1000):#compute autocorrelation of E(t)  ss1:=seq(x2(i),i=1..N):  ss2:=seq(E2(i),i=1..N):  ColumnGraph(Rx,'color'="Black",'style'='polygon');  ColumnGraph(RE,'color'="Red",'style'='polygon');  plot(<ss1>,style=point,symbol=point, color=blue);  plot(<ss2>,style=point,symbol=point,color=green);   Eav:= (Etot/N)*nm^2; # <E> average step (4)    E2av:= (E2tot/N)*nm^4; # <E2>  Xav:=(Xtot/(N))*nm^(2);#<X> average  X2av:=(X2tot/(N))*nm^(4);#<X^(2)> average   CV:= (E2av-Eav^2)/(kB*T^2)/nm^2;

Just one thing that I am not sure of is the conditional, does it satisfy my official problem in my OP?

Here it is:

"perfroming individual steps by selecting an attmept to move x=+-a ,with equal probability, and choose whether to make the actual step according to energy change expected in that step."

## I am attaching below my code that I chan...

I am attaching below my code that I changed a little bit from the code in the link in my OP, I added autocorrelation functions and the graphs of theirs and of the displacement and energy.

I get the following error:

Error, cannot determine if this expression is true or false: 0.2071500000e-4*x2^2 <= 0. or 149.5449674 < exp(-0.5000000000e-2*x2^2)

It seems that something is wrong in the if conditional, how to fix it?

Also, the values of <E^2>,<E>,<x>,<x^2>,C_V are all zero, which doesn't seem to me to be probable.

Can anyone help me with my code?

Thanks.

## @vv where can I find this Poisson f...

@vv where can I find this Poisson formula for 2d Poisson equation?

In Evans' book on PDE second edition I have the following solution for n=2:

u(x) = -1/2pi \int_{\mathbbb{R}^2} log(|x-y|)f(y)dy

## I need some assistance....

@tomleslie , hi tom; I changed a little bit your code (the natural logs I changed to log10, tval changed to 5, and I am now looking at the error norm between cos(x+t) and the solution of the PDE:

diff(v(x, t), t) = diff(v(x, t), x, x)-sin(x+t)+cos(x+t)

D[1](v)(0,t)=-sin(t),
D[1](v)(1,t)=-v(1, t)^4+(cos(1+t))^4-sin(1+t),
v(x,0)=cos(x)

the analytical solution of this PDE is cos(x+t))

Now in order to calculate the slope of the log10 log10 graph I use:

with(CurveFitting):
LeastSquares([[log10(1/16), log10(E[1])], [log10(1/32),log10(E[2])],[log10(1/64), log10(E[3])],[log10(1/128),log10(E[4])],[log10(1/256),log10(E[5])]], v);

I get a slope that is smaller than 2 while if I change the boundary conditions to:

D[1](v)(0,t)=-sin(t),
D[1](v)(1,t)=-sin(1+t),
v(x,0)=cos(x)

(its analytical solution is also cos(x+t))

I get a slope which is approximately 2.

My colleague told me that theoretically the slope should be 2, why is that if I change to the first boundary conditions that the slope of the log log graph decreases so drastically?

The analytical solution of the two problems is identical, isn't it?

pdeProb_(1).mw

P.S

Forgot to add that now we are comparing the numerical solutions of the PDEs with the analytical solution which in this case is cos(x+t).

## @tomleslie  Hi, if I would like to ...

Hi, if I would like to use leastsquares command and plot the function log ||E||_h as a function of log h, how would I change your code?

I mean we have:

pA[2]:= plot( [ seq( [ log(hvals[j]), log(E[j]) ],
j=1..numelems(hvals)
)
],
labels=["space\\timestep", "Norm"],
labelfont=[times, bold, 16],
labeldirections=[horizontal, vertical],
title="Norm vs step (log10-log10)",
titlefont=[times, bold, 20]
):

plots:-display(pA);

So I would like to generate a least squares approximation with the xdata being log h=log 1/16 ,..., log 1/256

and the ydata would be the corresponding value of log ||E||_h at the points of h=1/16 ,...,1/256

and then to plot the least square approximation.

So the least squares part would be:

LeastSquares([[log(1/16), log(E[1])], [log(1/32),log(E[2])],[log(1/64), log(E[3])],[log(1/128,log(E[4])],[log(1/256),log(E[5])]], v);

Now, how to plot this least squares function?

as log ||E||_h vs log h.

## @tomleslie thanks Tom. I appreciate...

@tomleslie thanks Tom.

I appreciate your help, if I were to use the same analysis just instead of the error between the two PDEs, only to calculate the norm of the solution of the first PDE, all I need to change in your code is this line:

E[p]:= sqrt
( (sol[1](j*h))^2,
j=0..1/h
)
):

Or do I need to change also other lines of your code?

## @tomleslie the correct form syntact...

@tomleslie the correct form syntactically is:

## @Markiyan Hirnyk yes I asked on PDE...

@Markiyan Hirnyk yes I asked on PDEs, I forgot that I should be using 2^(-k) here.

## @tomleslie thanks, so much syntax t...

@tomleslie thanks, so much syntax to learn to use here.

## @tomleslie I am sorry that I am ask...

@tomleslie I am sorry that I am asking another question.

But how can I plug a value into mf after the substituion of the numerical paramters, I mean I want to plug in the IBC after the PDEs values of u(x,0) = 1+mf(at t=0) and u(x,0) = 1+myPoly(at t=0) respectively, how to implement this?

I found how to do this for mf, just subs[eval](t=0,mf)=:mf1 I am sure it will work also for the polynomial.

## @tomleslie , I have another request...

@tomleslie , I have another request, how do I use the coefficients of the two approaches (of a*exp(-b*t)-c*t and of the sixth polynomial) in any other further commands; i.e I want to use mf as an explicit function with the numerical coefficients obtained before and the same with the polynomial in further operations.

Here's an example of what I want to do:

PDE2 := diff(u(x, t), t) = diff(u(x, t), x, x)+sin(x+t)-cos(x+t)-(diff(mf, t));

IBC2 := (D[1](v))(0, t) = 0, (D[1](v))(1, t) = -0.65e-4*(v(1, t)+mf)^4, v(x, 0) = 1+a

pds2 := pdsolve(PDE2, [IBC2], numeric, time = t, range = 0 .. 1, spacestep = 0.1e-1, timestep = 0.1e-1, errorest = true)*pds2:-plot(x = 1, t = 0 .. 1)

How to make the polynomial and a*exp(-b*t)-c*t with the numerical coefficients available to me?

## @tomleslie thanks. The polynomial s...

@tomleslie thanks.

The polynomial seems like a close to accurate approximation.

## @tomleslie are there any other meth...

@tomleslie are there any other methods? how about polynomial interpolation?

I don't understand how do I extract the coefficients a,b,c from the fitting a*exp(-b*t)-c*t?

## @Preben Alsholm Hi, I am reading th...

@Preben Alsholm Hi, I am reading the book: Finite elements using maple.

They give a code for collocation approximation to a simple ODE with weights, can we use this approach here with my nonlinear integral equation?

with 100 collocation points x[i]?

Here's a preview of the book:

on pages 79-80 there's the code.

A description of collocation method is on page 68.

## @tomleslie my programming skills ar...

@tomleslie my programming skills are weak, I know, I am trying to ammend it.

 1 2 3 4 Page 2 of 4
﻿