

In [ ]:

    
using BasisFunctions
using FrameFun
using Plots
using DomainSets



In [ ]:

    
FE = FrameFun
BA = BasisFunctions

1D Differential Equation (Poisson)



In [ ]:

    
B = FourierBasis(41,-1,1)
Dom = -0.5..0.5

Homogenous Dirichlet Boundary Condition

$$\begin{align}p''(x)&=f(x), &x\in\Omega\\ p(x)&= 0, &x \in \delta\Omega\end{align}$$

Boundary conditions are defined by an operator and a function. The solution F satisfies diff*F = df on the domain boundary.



In [ ]:

    
diff = IdentityOperator(B)
BC = DirichletBC(x->0);

Differential equations are defined by an operator, a function and a boundary condition. The solution F satisfies Diff*F=f in the interior of the domain.



In [ ]:

    
fD = x->x;
Diff = differentiation_operator(B)*differentiation_operator(B)
DE = DiffEquation(B,Dom,Diff,fD, (BC,BC));



In [ ]:

    
FD = solve(DE)



In [ ]:

    
# Exact solution
solD = x->x^3/6 - x/24;



In [ ]:

    
plot(FD,layout=4,title="Solution")
plot!(FD'',subplot=2,title="Second derivative")
plot!(FD,solD,subplot=3,title="Solution error")
plot!(FD'',fD,subplot=4,title="Derivative error")

Homogenous Neumann Boundary Condition

$$\begin{align}p''(x)&=f(x), &x\in\Omega\\ p'(x)&= 0, &x \in \delta\Omega\end{align}$$



In [ ]:

    
diff = differentiation_operator(B)
NeumannBC(x->0.)



In [ ]:

    
fN = x->x;
Diff = differentiation_operator(B)*differentiation_operator(B)
DE = DiffEquation(B,Dom,Diff,fN, (BC,));



In [ ]:

    
FN = solve(DE)



In [ ]:

    
# Exact solution
solN = x->x^3/6-x/8



In [ ]:

    
plot(FN,layout=4,title="Solution")
plot!(FN'',subplot=2,title="Second derivative")
plot!(FN,solN,subplot=3,title="Solution error")
plot!(FN'',fN,subplot=4,title="Derivative error")

2D Differential Equation

2D experiments may take some time

Dirichlet on Annulus (Laplace)

$$\begin{align} \Delta p(x,y)&=0 &(x,y)\in\Omega\\ p(x,y) &= df(x,y) &(x,y)\in\delta\Omega \end{align}$$



In [ ]:

    
B2 = FourierBasis(31,-1.0,1.0)⊗FourierBasis(31,-1.0,1.0)
D2 = disk(0.8)\disk(0.3)\cube([-0.15,-1.0],[0.15,0.0])



In [ ]:

    
diff2 = IdentityOperator(B2)
df2D = (x,y) -> x-y;
BC2 = DirichletBC(df2D,euclideanspace(Val{2}()))



In [ ]:

    
f2D = (x,y)->0;
Diff2 = differentiation_operator(B2,(2,0))+differentiation_operator(B2,(0,2))
DE2 = DiffEquation(B2,D2,Diff2,f2D, (BC2,));



In [ ]:

    
F2D = solve(DE2,solver=DirectSolver)



In [ ]:

    
heatmap(F2D)

Homogenous Neumann on semi-periodic strip (Poisson)

$$\begin{align} \Delta p(x,y)&=f(x,y) &(x,y)\in\Omega\\ \frac{\delta p}{\delta y}(x,y) &= 0 &(x,y)\in\delta\Omega \end{align}$$



In [ ]:

    
B2 = FourierBasis(21,-1.0,1.0)⊗FourierBasis(21,-1.0,1.0)
D2 = cube([-1.0,-0.5],[1.0,0.5])



In [ ]:

    
diff2 = differentiation_operator(B2,(0,1))
df2N = (x,y)->0;
BC2 = NeumannBC(df2N,euclideanspace(Val{2}()))



In [ ]:

    
f2N = (x,y)->cos(2*pi*(x+y));
Diff2 = differentiation_operator(B2,(2,0))+differentiation_operator(B2,(0,2))
DE2 = DiffEquation(B2,D2,Diff2,f2N, (BC2,));



In [ ]:

    
F2N = solve(DE2,solver=AZSolver)



In [ ]:

    
heatmap(F2N)

Helmholtz



In [ ]:

    
using StaticArrays
B2 = FourierBasis(41,-1.0,1.0)⊗FourierBasis(41,-1.0,1.0)
D2 = disk(0.75)\disk(0.2,SVector(0.3,-0.3))



In [ ]:

    
diff2 = IdentityOperator(B2)
df2H = (x,y)->0;
BC = NeumannBC(df2H,euclideanspace(Val{2}()))



In [ ]:

    
f2H = (x,y)->exp(-200*((x+0.3)^2+(y-0.3)^2))
Diff2 = differentiation_operator(B2,(2,0))+differentiation_operator(B2,(0,2))+1000*IdentityOperator(B2)
DE2 = DiffEquation(B2,D2,Diff2,f2H, (BC2,));



In [ ]:

    
F2H = FrameFun.solve(DE2,solver=AZSolver)



In [ ]:

    
heatmap(F2H)



In [ ]: