@vv thanks.

@vv How does it work for every t>0 ?

@vv Is there a way for maple to show me the calculation explicitly?

@vv Yes, I agree, I'll ask my adviser for further instructions.

@vv I agree, but in my case the argument depends on other terms that are explicit functions of t, so if we maximize the all whole argument, then if it were to depend on p_0(t) then the max would depend on max p_0(t).

Here are a few examples to what I mean:

\sup_{t\in [0,1/epsilon} |t+p_0(t)| = 1/\epsilon +\sup_{t\in [0,1/\epsilon)}|p_0(t)|

\sup_{t\in [0,1/epsilon) |t-p_0(t)| = 1/\epsilon +\max{-\inf_{t\in [0,1/epsilon)} |p_0(t)|,\sup_{t\in [0,1/epsilon) |p_0(t)|}

Both cases are legitimate instances of what I am searching, can it be done in maple 2017 with my case?

Obviously the dependence of \delta_1(\epsilon) will be on \sup_{t\in [0,1/\epsilon)} or \inf_{t\in [0,1/\epsilon} of |p_0(t)|

but I want to see this dependence explciitly.

@tomleslie Thanks for that.

I tried to compute \delta1(\epsilon) but I get that it doesn't compute the max with respect to t.

I think the x dependence gets cancelled.

I am attaching a file with what I tried to do here.

Perturbation.mw

Can you help me with defining the correct \delta(\epsilon) that it will compute it correctly?

Thanks in advance!

I think I understand how to define most of the definitions here, but I am not sure how to define that p_0(t) indeed will be taken as a function of t and the max function will regard it as such.

Here are the functions that I've defined so far:

f := proc (x, y1, y2, z) x+y1^2*z end proc

h := proc (x, y1, y2) y1^2 end proc

q_0 := proc (s) exp(-s)*(q0+int(h(x0, cos(r), sin(r))-p0, r = 0 .. s)) end proc

g := proc (x, t) Limit(int(f(x, cos(s), sin(s), p_0(t)+q_0(s)), s = 0 .. T)/T, T = infinity) end proc

Now I need to compute the maximum such that p_0(t) will be regarded as an unkown function of t.

BTW, when I try to compute g(x,t) I get the following answer:

              /    1    //           2                           
     lim      |- ------ \\1360 cos(T)  q0 - 2720 cos(T) sin(T) q0
T -> infinity \  6800 T                                          

   - 680 cos(2 T) T p0 + 1360 sin(2 T) T p0

   - 3400 cos(T) sin(T) exp(T) p_0(t) + 340 cos(2 T) T

   + 408 cos(2 T) p0 - 680 sin(2 T) T + 544 sin(2 T) p0

   + 100 cos(4 T) + 25 sin(4 T) + 136 cos(2 T) - 102 sin(2 T)

   - 3400 exp(T) p_0(t) T - 6800 x T exp(T) + 2992 p0 exp(T)

   - 4080 q0 exp(T) - 3400 T p0 - 1936 exp(T) + 1700 T - 3400 p0

                   \        \\
   + 2720 q0 + 1700/ exp(-T)/|
                             /

Why doesn't it compute the limit?

THanks.

@Rouben Rostamian  yes you are right, as I formed it doesn't seem meaningful.

I'll repost my question after I refine exactly what I am looking for.

@acer can you give me exact instructions how to fix this? I mean what exactly did you do to fix this in your system?

Thanks!

So I think I succesully installed both maple 2017.3 and the toolbox for matlab.

The problem is while running maple through root, i.e running the root in the terminal, I still get the following error:

Exception in thread "Request id 1" java.lang.InternalError: XXX0 profile[1]: GL3bc -> profileImpl GL4bc !!! not mapped
    at com.jogamp.opengl.GLProfile.computeProfileMap(GLProfile.java:2071)
    at com.jogamp.opengl.GLProfile.initProfilesForDeviceCritical(GLProfile.java:1924)
    at com.jogamp.opengl.GLProfile.initProfilesForDevice(GLProfile.java:1875)
    at com.jogamp.opengl.GLProfile.initProfilesForDefaultDevices(GLProfile.java:1843)
    at com.jogamp.opengl.GLProfile.access$000(GLProfile.java:80)
    at com.jogamp.opengl.GLProfile$1.run(GLProfile.java:230)
    at java.security.AccessController.doPrivileged(Native Method)
    at com.jogamp.opengl.GLProfile.initSingleton(GLProfile.java:216)
    at com.jogamp.opengl.GLProfile.getProfileMap(GLProfile.java:2297)
    at com.jogamp.opengl.GLProfile.get(GLProfile.java:988)
    at com.jogamp.opengl.GLProfile.get(GLProfile.java:1004)
    at com.maplesoft.mathdoc.view.plot.Plot3DCanvasView.<clinit>(Unknown Source)
    at com.maplesoft.worksheet.view.WmiWorksheetView.readProperties(Unknown Source)
    at com.maplesoft.worksheet.view.WmiWorksheetView.wmiWorksheetView_init(Unknown Source)
    at com.maplesoft.worksheet.view.WmiWorksheetView.<init>(Unknown Source)
    at com.maplesoft.worksheet.view.WmiWorksheetView.<init>(Unknown Source)
    at com.maplesoft.worksheet.application.WmiWorksheetWindowManager.createWorksheet(Unknown Source)
    at com.maplesoft.worksheet.application.WmiWorksheetWindowManager.createWorksheet(Unknown Source)
    at com.maplesoft.worksheet.application.WmiWorksheetWindowManager.createWorksheet(Unknown Source)
    at com.maplesoft.worksheet.application.WmiWorksheetWindowManager.createWorksheet(Unknown Source)
    at com.maplesoft.worksheet.application.WmiWorksheetStartupStrategy.createFirstWorksheet(Unknown Source)
    at com.maplesoft.worksheet.application.WmiWorksheetStartupStrategy.openFirstWorksheet(Unknown Source)
    at com.maplesoft.worksheet.application.WmiGenericStartupStrategy.doStartup(Unknown Source)
    at com.maplesoft.worksheet.application.WmiWorksheetStartupStrategy.doStartup(Unknown Source)
    at com.maplesoft.application.Maple.doStartup(Unknown Source)
    at com.maplesoft.application.Application.startup(Unknown Source)
    at com.maplesoft.application.ServerProtocol$StartApplicationHandler.processCommand(Unknown Source)
    at com.maplesoft.application.ServerProtocol.executeCommand(Unknown Source)
    at com.maplesoft.application.ServerProtocol.processNextStep(Unknown Source)
    at com.maplesoft.application.ExchangeProtocol.executeProtocol(Unknown Source)
    at com.maplesoft.application.ApplicationManager$Listener.run(Unknown Source)
    at java.lang.Thread.run(Thread.java:745)
And maple initial app window gets stuck, I searched google for the first line error:

Exception in thread "Request id 1" java.lang.InternalError: XXX0 profile[1]: GL3bc -> profileImpl GL4bc !!! not mapped

But I am not sure that's relevant for my problem in maple.

Can anyone help me with this?

@Carl Love 

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

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

Thanks in advance, I appreciate your help.

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