```julia
using MDBM
using PyPlot;
pygui(true);
```

# 1. Implcit equation with constraint
\begin{align}
 x,y & \in \left[ -5,5 \right] \\
x^2+y^2-2^2 & =0 \\ 
x-y & >0
\end{align}


```julia
function foo(x,y)
    x^2.0+y^2.0-4.0^2.0
end
function c(x,y)
    x-y
end

ax1=Axis([-5,-2.5,0,2.5,5],"x")
ax2=Axis(-5:2:5.0,"b")

mymdbm=MDBM_Problem(foo,[ax1,ax2],constraint=c)
iteration=6 #number of refinements (resolution doubling)
solve!(mymdbm,iteration)

#evaluated points
x_eval,y_eval=getevaluatedpoints(mymdbm)

#solution points
x_sol,y_sol=getinterpolatedsolution(mymdbm)

fig = figure(1);clf()
scatter(x_eval,y_eval,s=2)
scatter(x_sol,y_sol,s=4);

# plot the constraint
mymdbm_c=MDBM_Problem(c,[ax1,ax2])
solve!(mymdbm_c,iteration)
x_sol,y_sol=getinterpolatedsolution(mymdbm_c)
scatter(x_sol,y_sol,s=4);
```


```julia
fig = figure(101);clf()
myDT1=connect(mymdbm);
for i in 1:length(myDT1)
    dt=myDT1[i]
    P1=getinterpolatedsolution(mymdbm.ncubes[dt[1]],mymdbm)
    P2=getinterpolatedsolution(mymdbm.ncubes[dt[2]],mymdbm)
    plot([P1[1],P2[1]],[P1[2],P2[2]], color="k")
end 
```

# 2. System of implicit equations
## parameter space

\begin{equation}
x,y,z \in \left[ -2,2 \right] 
\end{equation}

## Sphere 1
Single implicit eqution
\begin{equation}
x^2+y^2+z^2-1 =0 
\end{equation}

## Sphere 2
Single implicit eqution
\begin{equation}
(x-0.5)^2+(y-0.5)^2+(z-0.5)^2-1 =0 
\end{equation}

## Intersection of two spheres
System of implicit equtions
\begin{align}
x^2+y^2+z^2-1 & =0  \\
(x-0.5)^2+(y-0.5)^2+(z-0.5)^2-1 & =0 
\end{align}


```julia
using LinearAlgebra
axes=[-2:2,-2:2,-2:2]
fig = figure(2);clf()

#Sphere1
fS1(x...)=norm([x...],2.0)-1
Sphere1mdbm=MDBM_Problem(fS1,axes)
solve!(Sphere1mdbm,4)
a_sol,b_sol,c_sol=getinterpolatedsolution(Sphere1mdbm)
plot3D(a_sol,b_sol,c_sol,linestyle="", marker=".", markersize=1);

#Sphere2
fS2(x...)=norm([x...] .- 0.5,2.0) -1.0
Sphere2mdbm=MDBM_Problem(fS2,axes)
solve!(Sphere2mdbm,4)
a_sol,b_sol,c_sol=getinterpolatedsolution(Sphere2mdbm)
plot3D(a_sol,b_sol,c_sol,linestyle="", marker=".", markersize=1);

#Intersection
# fS12(x...)=[fS1(x...),fS2(x...)]

function fS12(x,y,z)
    fS1(x,y,z),fS2(x,y,z)
end
Intersectmdbm=MDBM_Problem(fS12,axes)
solve!(Intersectmdbm,6)
a_sol,b_sol,c_sol=getinterpolatedsolution(Intersectmdbm)
plot3D(a_sol,b_sol,c_sol,color="k",linestyle="", marker=".", markersize=2);
```

# 3. Non-smooth problem
## Mandelbrot set

Convergence test of series $z_{i+1}=z_i^2+c$ for $z_0=0$, $c \in  \mathbb{C}$

In the evaluation $c=x + i y$

If $\| z_i \| >2$, then the series will diverge.



```julia
function mandelbrot(x,y)    
    c=x+y*1im
    z=Complex(0)
    k=0
    maxiteration=1000
    while (k<maxiteration && abs(z)<4.0)
            z=z^2+c
            k=k+1
    end
    
    abs(z)-2.0
end

Mandelbrotmdbm=MDBM_Problem(mandelbrot,[-5:2,-2:2])
solve!(Mandelbrotmdbm,8)
real_c_sol,imag_c_sol=getinterpolatedsolution(Mandelbrotmdbm)
fig = figure(3);clf()
plot(real_c_sol,imag_c_sol,linestyle="", marker=".", markersize=1);
```

# 4. Continuation behavioue
## 4.1 Problem:
\begin{equation}
 x,y \in \left[ -2,2 \right]
\end{equation}

\begin{equation}
\| \mathbf{x} +1.5\|_{1.7}+\sin(5x_1)-2 =0 
\end{equation}


```julia
using LinearAlgebra
mymdbm=MDBM_Problem((x...)->norm([x...].+ 1.5,1.7)+sin(x[1]*5)-2.0,[-2:2,-2:2])

interpolate!(mymdbm,interpolationorder=1)
for k=1:5
    refine!(mymdbm)
    interpolate!(mymdbm,interpolationorder=1)
end
checkneighbour!(mymdbm)
interpolate!(mymdbm,interpolationorder=1)

#evaluated points
a_eval,b_eval=getevaluatedpoints(mymdbm)

a_sol,b_sol=getinterpolatedsolution(mymdbm)

fig = figure(4)
plot(a_eval,b_eval,linestyle="",  marker=".",markersize=1)
plot(a_sol,b_sol,linestyle="",  marker=".",markersize=1);
```

## 4.2 Exploring the missing component

Due to the coarse initial mesh, some part of the solution is missing.


A continuation-like exploration of the missing component can be performed by checking the neighbouring n-cubes.

It is also clear, that the range of the initial mesh do not cover the object. The actual grid can be prepend and append also.

Note, that the extended grid is used only for the continuation of the detected segments! The initialized new grid points are not evaluated.<br>
In this example there is a closed curve of soultion around $x_1=-4.1$ which is lost!


```julia
# # extension with the same resolution:
# axesextend!(mymdbm,1,prepend=mymdbm.axes[1].ticks[1:end-1] .+(mymdbm.axes[1].ticks[1]- mymdbm.axes[1].ticks[end]));
# axesextend!(mymdbm,2,prepend=mymdbm.axes[2].ticks[1:end-1] .+(mymdbm.axes[2].ticks[1]- mymdbm.axes[2].ticks[end]));

#extension with different reselution
axesextend!(mymdbm,1,prepend=-6.2:0.1:-2.2);
axesextend!(mymdbm,2,prepend=-6.2:0.05:-2.2,append=2.2:0.2:3);


checkneighbour!(mymdbm)
interpolate!(mymdbm,interpolationorder=1)

#evaluated points
a_eval,b_eval=getevaluatedpoints(mymdbm)

#solution points
a_sol,b_sol=getinterpolatedsolution(mymdbm)


# scatter(a_eval,b_eval,markersize=1)
# scatter!(a_sol,b_sol,size = (800, 800),markersize=2,
#      xticks = mymdbm.axes[1].ticks , yticks = mymdbm.axes[2].ticks, gridalpha=0.8)

fig = figure(4)
clf()
plot(a_eval,b_eval,linestyle="",  marker=".",markersize=1)
plot(a_sol,b_sol,linestyle="",  marker=".",markersize=1);
```

## 4.3 Further refinement
The prepended and appended grid has a poor resolution, thus further refinement can be used to increase the resolution.


```julia
for k=1:2
    refine!(mymdbm)
    interpolate!(mymdbm,interpolationorder=1)
end
checkneighbour!(mymdbm)
interpolate!(mymdbm,interpolationorder=1)

#evaluated points
a_eval,b_eval=getevaluatedpoints(mymdbm)

#solution points
a_sol,b_sol=getinterpolatedsolution(mymdbm)


# scatter(a_eval,b_eval,markersize=1)
# scatter!(a_sol,b_sol,size = (800, 800),markersize=2,
#  xticks = mymdbm.axes[1].ticks , yticks = mymdbm.axes[2].ticks, gridalpha=0.8)

fig = figure(4)
clf()
plot(a_eval,b_eval,linestyle="",  marker=".",markersize=1)
plot(a_sol,b_sol,linestyle="",  marker=".",markersize=1);
```

# 5. Constraint only
If only a constraint if provided, then a dens point cloud will be provided wich could be used for later computataions.
## Problem:
\begin{equation}
 x,y \in \left[ -5,2 \right]
\end{equation}

\begin{equation}
\| \mathbf{x} +1.5\|_{1.7}+\sin(5x_1)-2 <0 
\end{equation}

Note, in this case the functinos will be evaluated in all the gridpoint of the interior.


```julia
using LinearAlgebra

ax1=Axis(-5:3.0,"a")
ax2=Axis(-5:3.0,"b")

mymdbm=MDBM_Problem((x...)->[],[ax1,ax2],constraint=(x...) -> -(norm([x...].+ 1.5,1.7)+sin(x[1]*5)-2.0))
# mymdbm=MDBM_Problem((x...) -> -maximum([0.0,-(norm([x...].+ 1.5,1.7)+sin(x[1]*5)-2.0)]),[ax1,ax2])
interpolate!(mymdbm,interpolationorder=1)
for k=1:4
    refine!(mymdbm)
    interpolate!(mymdbm,interpolationorder=1)
end
checkneighbour!(mymdbm)
interpolate!(mymdbm,interpolationorder=1)

#solution points
a_sol,b_sol=getinterpolatedsolution(mymdbm)

F_sol=map((x,y)->x*x-y*y,a_sol,b_sol)
# scatter(a_sol,b_sol,F_sol,size = (500, 500))


fig = figure(5)
plot3D(a_sol,b_sol,F_sol,linestyle="",  marker=".",markersize=4);
```

# 6. Connection of point could
The generated points coloud can be connected by generatigh the eges between the neighbouring n-cubes.
## Problem:

\begin{align}
 x,y,z & \in \left[ -5,5 \right] \\
x^2+y^2+z^2-4 & =0 \\ 
z+y-\cos(3x) & >0
\end{align}


```julia
ax1=Axis(-3.0:3.0,"x")
ax2=Axis(-3.0:3.0,"y")
ax3=Axis(-3.0:3.0,"z")

function foo(x,y,z)
    [x^2+y^2+z^2-4]
end

function cont(x,y,z)
    [z+y-cos(3*x)-1]
end

mymdbm=MDBM_Problem(foo,[ax1,ax2,ax3],constraint=cont)
solve!(mymdbm,2)
        
#solution points
x_sol,y_sol,z_sol=getinterpolatedsolution(mymdbm);
fig = figure(10),clf()
plot3D(x_sol,y_sol,z_sol,linestyle="",  marker=".",markersize=4);

#connecting the neighbouring n-cubes
myDT1=connect(mymdbm);

#plot solution lines one-by-one
fig = figure(11);clf()
for i in 1:length(myDT1)
    dt=myDT1[i]
    P1=getinterpolatedsolution(mymdbm.ncubes[dt[1]],mymdbm)
    P2=getinterpolatedsolution(mymdbm.ncubes[dt[2]],mymdbm)
    plot3D([P1[1],P2[1]],[P1[2],P2[2]],[P1[3],P2[3]], color="k")
end 
```

# 7. Complex problem: 4 parameter, 2 equation, 1 constraint
## Problem:

\begin{align}
 x,y,z \in \left[ -3,3 \right]  & \quad r \in \left[ 1,2 \right] \\
x^2+y^2+z^2-r^2 & =0 \\ 
z+y-\cos(3x) & =0 \\
z-\sin(5y) & >0
\end{align}


```julia
ax1=Axis(-3.0:3.0,"x")
ax2=Axis(-3.0:3.0,"y")
ax3=Axis(-3.0:3.0,"z")
ax4=Axis(1.0:0.5:2.0,"r")

function foo(x,y,z,r)
    [x^2+y^2+z^2-r^2,
        z+y-cos(3*x)]
end

function cont(x,y,z,r)
   z-sin(5*y)
end

mymdbm=MDBM_Problem(foo,[ax1,ax2,ax3,ax4],constraint=cont)
solve!(mymdbm,3)

#solution points
x_sol,y_sol,z_sol,r_sol=getinterpolatedsolution(mymdbm);
fig = figure(20);clf()
plot3D(x_sol,y_sol,z_sol,linestyle="",  marker=".",markersize=4);
```

6226462123e92349566368fe926f6a706830542e
ipynb
Jupyter Notebook
examples/test.ipynb
arturgower/MDBM.jl
8c60ff4abffb308d5b1b8394d37ff8daa8640b53
MIT
2018-11-22T16:08:40.000Z
2021-12-07T18:29:07.000Z
examples/test.ipynb
arturgower/MDBM.jl
8c60ff4abffb308d5b1b8394d37ff8daa8640b53
MIT
2019-02-20T15:48:21.000Z
2021-12-15T20:19:02.000Z
examples/test.ipynb
arturgower/MDBM.jl
8c60ff4abffb308d5b1b8394d37ff8daa8640b53
MIT
2018-08-14T13:56:12.000Z
2021-12-09T13:43:35.000Z
Qwen/Qwen-72B
1. YES
2. YES