maroba / findiff

Python package for numerical derivatives and partial differential equations in any number of dimensions.
MIT License
420 stars 60 forks source link

Two boundary conditions on the same boundary #22

Closed PetterTh closed 1 year ago

PetterTh commented 4 years ago

Hi

How would you go about to set two boundary conditions on the same border? For instance in you example with the 1D forced harmonic oscillator with friction to have both u(0)=0 and du/dx (0) = 0.

Petter

maroba commented 4 years ago

Hi, what you describe would be an initial condition problem. Currently findiff only handles boundary value problems. For a given boundary point, you can only specify one condition, otherwise the system would be overdetermined.

PetterTh commented 4 years ago

Ok, thank you for the clarification. But it's not only an issue for initial value problems. Solving the equation for the Euler Beam (https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory#Static_beam_equation) require 4 boundary conditions in total with only two boundaries. Do you think it is something which will be possible to do in the future with this package?

maroba commented 4 years ago

Sorry for the late response. I haven't had time for this project lately. Regarding your question, I must say that the PDE functionality is actually a side feature of this package which actually is only about numerical derivatives. If I have some time, I may go forward into this direction, but please don't count on this in the near future.

dangeo84 commented 11 months ago

I believe for cases like solving the beam equation, you create ghost nodes that exist outside of the physical beam and apply the extra boundary conditions to those.

#### SIMPLY SUPPORTED BEAM ####
## Set up the grid
l = 10.0  ## beam length
dx = 0.1
n = int(l//dx + 2)
shape = (n, )
x = np.linspace(-dx, l+dx, shape[0])

## define variables
d = 0.75  ## circular section diameter
q = 10  ## Uniformly distributed load
EpIp = 32e6*np.pi*d**4/64.0 ## Flexural stiffness for 750mm circular element

## Create the PDE
L = EpIp*FinDiff(0, dx, 4)
f = -q*np.ones(shape)

## Apply the boundary conditions for a simply supported beam
bc = BoundaryConditions(shape)
bc[0] = (FinDiff(0,dx,2), 0) ## Enforce zero curvature at left-most ghost point
bc[1] = 0   ## Enforce zero displacement at left support

bc[-1] = (FinDiff(0,dx,2), 0)   ## As above, but for right-most ghost point
bc[-2] = 0  ## As above but for right support

## Solve
pde = PDE(L, f, bc)
y = pde.solve()

## Plot
plt.plot(x,y)
plt.gca().set_xlim(left=0,right=10)
plt.grid()