#################################################################################################### # This is a program to test whether a function satisfies a fully nonlinear parabolic equation. # We also have a slightly more efficient version. This one is easier to use. # # In order to run this, you need to have julia installed, and know how to use an interactive session # See https://julialang.org/ # # It requires the packages LinearAlgebra and SymEngine, as well as Distributed for a parallel computation. # See https://github.com/symengine/SymEngine.jl for instructions on how to install SymEngine # #################################################################################################### # In order to use this, open an interactive julia session and type # include("verifyD_simple.jl") # # Then, run these commands # u = u0() # sample_and_walk_function(u,100,2000,9.,12.) # # Here, u is the candiate solution. # 1000 is the number of sample initial points for the gradient flow. # 2000 is the number of iterations for the stochastic gradient flow for each initial point. # If the ratio is ever larger than 9., it performs 10x2000 more iterations. # It the ratio is ever larger than 12., it performs 100x2000 more iterations. # It returns the largest ratio is finds. # Note that the stochastic random walk stops if the ratio is ever larger than 1000. # Modify the parameters as desired # # We can also verify the singular solution to an elliptic PDE by Nadirashvili-Tkachev-Vladut with # the following code: # u = NVT() / sqrt(n2) # sample_and_walk_function(u,50,10000) # Interestingly, it returns a maximum ratio of approximately 9. See Remark 5.3 in the paper # # In order to verify the example by Charles Smart, use the following commands # change_dimension!(9) # u = (x[1]*x[5]*x[9]+x[2]*x[6]*x[7]+x[3]*x[4]*x[8]-x[1]*x[6]*x[8]-x[2]*x[4]*x[9]-x[3]*x[5]*x[7])/sqrt(n2) # sample_and_walk_function(u,50,10000) # #################################################################################################### # # This code is trivially parallelizable by just performing independent samples. # If you want to run a parallel computation, just uncomment the following 4 lines as well # as the code at the end of the file ###################BEGINNING OF PARALLEL CODE (UNCOMMENT IF DESIRED)############################### # using Distributed # # Set CPUs to be equal to the number of cores available in your computer. # # The following lines are unnecessary if you start julia with the parameter -p CORES # CPUs = 4 # addprocs(CPUs) # @everywhere begin #################END OF CODE FOR PARALLEL COMPUTATION############################################## dimension = 5 const tolerance = 1e-10 # below this we consider zero using SymEngine using LinearAlgebra # We define the symbols t, x1, x2, x3, ... t = symbols("t") x = Array{Basic,1}(undef,dimension) for i in 1:dimension x[i] = symbols("x$i") end n2 = reduce(+,map(x->x^2,x)) function NVT() return (x[1]^3 + 3x[1]/2 * (x[3]^2+x[4]^2-2*x[5]^2-2x[2]^2) + 3*sqrt(3)/2 * (x[2]*x[3]^2 - x[2]*x[4]^2 + 2*x[3]*x[4]*x[5])) end """ This function returns our current best candidate for a singular solution to a parabolic equation. This is the singular solution that we report in the paper. The largest ratio reported by sample_and_walk_function so far with this function is 13.667 """ function u0() u = NVT()/ (sqrt(n2+t^2)-t) + (1/12) * NVT() return u end """ Use this function if you ever want to test for singular solutions in any dimension other than 5 """ function change_dimension!(d::Int) global dimension = d global x = Array{Basic,1}(undef,dimension) for i in 1:dimension x[i] = symbols("x$i") end global n2 = reduce(+,map(x->x^2,x)) end """ This function tells us the ratio between the positive and negative parts of (vt,DH) """ function pucci_ratio(vt::Real,DH::Array{<:Real,2}) e = eigvals(DH) Pp = 0. Pm = 0. if vt>0 Pp += vt else Pm -= vt end for ev in e try if ev>0 Pm += ev else Pp -= ev end catch println("Non symmetric matrix :",DH-transpose(DH)) return pucci_ratio(vt,(DH+transpose(DH))/2) end end if (Pp < tolerance)&&(Pm < tolerance) return 1 end return max(Pp/Pm, Pm/Pp) end """ This is the function that performs the stochastic gradient flow. """ function pair_walk_sample(lut::Function, lD2u::Array{<:Function,2}, st1::Real, sx1::Array{<:Real,1}, st2::Real, sx2::Array{<:Real,1}, samples::Int=100, step_size::Real=0.004, noise::Real=0.) m = 0 notmoving = 0 for i in 1:samples notmoving += 1 if (noise>0. && rand()1) nt1, nx1, nt2, nx2 = shake_pair2(st1, sx1, st2, sx2, step_size/notmoving) else nt1, nx1, nt2, nx2 = st1, sx1, st2, sx2 end ut1 = lut(nt1,nx1...) ut2 = lut(nt2,nx2...) D2u1 = zeros(dimension,dimension) D2u2 = zeros(dimension,dimension) for i in 1:dimension for j in 1:dimension D2u1[i,j] = lD2u[i,j](nt1,nx1...) D2u2[i,j] = lD2u[i,j](nt2,nx2...) end end nm = pucci_ratio(ut1-ut2,D2u1-D2u2) if (nm>m) m = nm st1, sx1, st2, sx2 = nt1, nx1, nt2, nx2 # println(m, " \tSteps: ", notmoving," \t step size:", step_size) #if (m>13) # println("t1= ",st1, " x1=", sx1) # println("t2= ",st2, " x2=", sx2) #end step_size = min(max(500*step_size / notmoving, 0.00001),0.5) notmoving = 0 #if (nm > 12.4) # println("Best point: ",st1, ", ", sx1, ", ", st2, ", ", sx2) #end #if (nm > bm) # global bm = nm # global bx1 = nx1 # global bt1 = nt1 # global bx2 = nx2 # global bt2 = nt2 #end if (m>1000) break end end end return m, st1, sx1, st2, sx2 end function shake_pair(t1::Real, x1::Array{<:Real,1}, t2::Real, x2::Array{<:Real,1}, step_size::Real) mt = (t1+t2)/2 mx = (x1+x2)/2 dt = (t1-t2)/2 dx = (x1-x2)/2 nt = mt + step_size * (rand()-0.5) if (nt>0) nt=0. end if (nt<-1) nt=-1. end nx = zeros(dimension) for i in 1:dimension nx[i] = mx[i] + step_size * (rand()-0.5) end rd = sqrt(reduce(+,map(x->x^2,dx))+dt^2) # nr = max(rd + step_size * (rand()-0.5) / 10,0.0001) nr = max(rd * (1 + 20*step_size * (rand()-0.5)),0.00001) st = dt/rd + step_size * (rand()-0.5) sx = zeros(dimension) for i in 1:dimension sx[i] = dx[i]/rd + step_size * (rand()-0.5) end nt1 = nt + nr*st nx1 = nx + nr*sx nt2 = nt - nr*st nx2 = nx - nr*sx nt1 = max(-1.,min(nt1,0.)) nt2 = max(-1.,min(nt2,0.)) normx1 = reduce(+,map(x->x^2,nx1)) if (normx1>1) r = sqrt(normx1) nx1 /= r end normx2 = reduce(+,map(x->x^2,nx2)) if (normx2>1) r = sqrt(normx2) nx2 /= r end return nt1, nx1, nt2, nx2 end """ This is a function that, given a par of points t1,x1 and t2,x2, returns a random pair of points nearby. It is used as the random step for the stochastic gradient flow. It moves the middle point some distance comparable to step_size. The difference between them is moved a distance comparable to step_size times the current distance. This is a way to handle the discontinuity of the ratio when t1,x1 == t2,x2. """ function shake_pair2(t1::Real, x1::Array{<:Real,1}, t2::Real, x2::Array{<:Real,1}, step_size::Real) mt = (t1+t2)/2 mx = (x1+x2)/2 dt = (t1-t2)/2 dx = (x1-x2)/2 nt = mt + step_size * (rand()-0.5) if (nt>0) nt=0. end if (nt<-1) nt=-1. end nx = zeros(dimension) for i in 1:dimension nx[i] = mx[i] + step_size * (rand()-0.5) end rd = sqrt(reduce(+,map(x->x^2,dx))+dt^2) # We do not want the two points of the pair to be identical if rd < 0.00001 dt *= 2 dx *= 2 end factor = 1.5 ndt = dt + factor * rd * step_size * (rand()-0.5) ndx = zeros(dimension) for i in 1:dimension ndx[i] = dx[i] + factor * rd * step_size * (rand()-0.5) end nt1 = nt + ndt nx1 = nx + ndx nt2 = nt - ndt nx2 = nx - ndx nt1 = max(-1.,min(nt1,0.)) nt2 = max(-1.,min(nt2,0.)) normx1 = reduce(+,map(x->x^2,nx1)) if (normx1>1) r = sqrt(normx1) nx1 /= r end normx2 = reduce(+,map(x->x^2,nx2)) if (normx2>1) r = sqrt(normx2) nx2 /= r end return nt1, nx1, nt2, nx2 end """ This is the main function to check whether a candidate function u solves some fully nonlinear parabolic equation. example usage: u = u0() sample_and_walk_function(u,1000,2000,9.,12.) Here, u is the candiate solution. 1000 is the number of sample initial points for the gradient flow. 2000 is the number of iterations for the stochastic gradient flow for each initial point. If the ratio is ever larger than 9., it performs 10x2000 more iterations. It the ratio is ever larger than 12., it performs 100x2000 more iterations. It returns the largest ratio is finds. Note that the stochastic random walk stops if the ratio is ever larger than 1000. """ function sample_and_walk_function(u::Basic, samples::Int, iterations::Int, threshold1=Inf, threshold2=Inf) ut = diff(u,t) D2u = [diff(u,x[a],x[b]) for a in 1:dimension, b in 1:dimension] lut = lambdify(ut,vcat([t],x)) lD2u = lambdify.(D2u,Ref(vcat([t],x))) maxr = 0 for i in 1:samples t1 = -rand() t2 = -rand() x1 = zeros(dimension) x2 = zeros(dimension) for j in 1:dimension x1[j] = 2*rand()-1 x2[j] = 2*rand()-1 end norm1 = sqrt(reduce(+,map(x->x^2,x1))) norm2 = sqrt(reduce(+,map(x->x^2,x2))) nn1 = rand()+10*tolerance nn2 = rand()+10*tolerance x1 *= nn1/norm1 x2 *= nn2/norm2 r, st1, sx1, st2, sx2 = pair_walk_sample(lut,lD2u,t1,x1,t2,x2,iterations) if r>threshold1 print("Iterating 10x steps. ") r, st1, sx1, st2, sx2 = pair_walk_sample(lut,lD2u, st1, sx1, st2, sx2,10*iterations) if r>threshold2 println("iterating 100x steps. ") r, st1, sx1, st2, sx2 = pair_walk_sample(lut,lD2u, st1, sx1, st2, sx2,100*iterations) println("point1: ", round(st1,digits=3),", ",round.(sx1,digits=3)) println("point2: ", round(st2,digits=3),", ",round.(sx2,digits=3)) print("------->") end end println("Max so far ",maxr,". This sample: ",r) if (r>maxr) maxr = r println("*******************************************^^^^^^^^^^^^^") end end return maxr end """ This is a utility function that given a function u and a point t0,x0, it returns the time derivative and the Hessian at that point. It is used to verify by hand the values of the ratio at the resulting sample points. """ function utD2u(u::Basic, t0::Real, x0::Array{<:Real,1}) ut = diff(u,t) D2u = [diff(u,x[a],x[b]) for a in 1:dimension, b in 1:dimension] lut = lambdify(ut,vcat([t],x)) lD2u = lambdify.(D2u,Ref(vcat([t],x))) ut = lut(t0,x0...) D2u = zeros(dimension,dimension) for i in 1:dimension for j in 1:dimension D2u[i,j] = lD2u[i,j](t0,x0...) end end return ut, D2u end ########################################################################################################## #### Uncomment the following code to run a parallel computation. #### Once this code is uncommented, an interactive julia session is no longer required. The computation #### can be carried out from command line simply by executing > julia verifyD_simple.jl #### It should return a maximum ratio between 13.4 and 13.7 at the end. ########################################################################################################## # end #@everywhere begin # # total_samples = 1000 # iterations = 2500 # println("Starting computation") # max_ratio = @distributed (max) for i in 1:CPUs # u = u0() # 1+1 # sample_and_walk_function(u,Int(round(total_samples/CPUs)),iterations,9.,12.) # end # println("The maximum ratio I found was ",max_ratio) ##########################################################################################################