# poisson.jl -- Solve the Poisson Equation # # del^2 y = f for x such that domain(x)>c # y = psi for x such that domain(x)=c using SparseArrays, LinearAlgebra, Plots xmin=[-2,-2] xmax=[2,2] mres=100 # specify the boundary condition function psi(x) return x[1]^3+exp(x[1])*x[1]*cos(x[2])+5 end # specify the distribution function f(x) return 6*x[1]+2*exp(x[1])*cos(x[2]) end # create spatial grid to fit the desired resolution if xmax[1]-xmin[1]1.0 count+=1 image[i1,i2]=count end end end return count end # create A and F using the 4-point stencil function stencil(count) xs=[]; ys=[]; as=[] fv=Vector{Float64}(undef,count) for i1=1:size(image)[1] for i2=1:size(image)[2] function look(j1,j2) if image[j1,j2]>0 push!(xs,image[i1,i2]) push!(ys,image[j1,j2]) push!(as,1.0) else fv[image[i1,i2]]-=psi(xi([j1,j2])) end end if image[i1,i2]>0 fv[image[i1,i2]]=dx^2*f(xi([i1,i2])) push!(xs,image[i1,i2]) push!(ys,image[i1,i2]) push!(as,-4.0) look(i1,i2-1) look(i1-1,i2) look(i1+1,i2) look(i1,i2+1) end end end return (sparse(Int.(xs),Int.(ys),Float64.(as)),fv) end # by construction the exact solution is psi function exact(count) xv=zeros(count) for i1=1:size(image)[1] for i2=1:size(image)[2] if image[i1,i2]>0 xv[image[i1,i2]]=psi(xi([i1,i2])) end end end return xv end # convert a solution back to the grid function physical(x) black=minimum(x) p=zeros(nu...) for i1=1:size(image)[1] for i2=1:size(image)[2] if image[i1,i2]>0 p[i1,i2]=x[image[i1,i2]] else p[i1,i2]=black end end end return p end # various plotting functions function pltdomain() heatmap(x1s,x2s,(x->min(1,x)).(image)',c=:hot) end function plterror() heatmap(x1s,x2s,physical(AP-EX)'*1000,c=:hot) end function pltexact() heatmap(x1s,x2s,physical(EX)',c=:hot) end function pltapprox() heatmap(x1s,x2s,physical(AP)',c=:hot) end # solve the problem and print the error count=enumerate() EX=exact(count) A,F=stencil(count) AP=A\F println(mres," ",norm(AP-EX)*dx)