Changing plot axis in heat-like equation

variational formulation
1

I’m solving a PDE with Dirac delta initial condition

\frac{\partial P}{\partial L}(L,\eta) = \frac{1}{L_{loc}}\frac{\partial}{\partial\eta}\bigg((\eta^2-1)\frac{\partial p}{\partial\eta}(L,\eta)\bigg), \quad \eta >1,

with initial condition P(L=0) = \delta(\eta-1). (\eta has been replaced with x). I discretize the length variable L via

\bigg(\frac{\partial P}{\partial L}\bigg)^{n+1} \approx \frac{P^{n+1}-P^n}{\Delta L},

the code below solves this problem as expected

"""
1d Fokker-Planck equation, spatially dependant diffusion process with
Dirac delta initial condition, Dirichlet and Neumann.
"""

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt

L = 2          # final L
num_steps = 10   # number of length steps
dl = L / num_steps # length step size
theta = 1 #Lloc constant term


# Create mesh and define function space
nx = 100
mesh = IntervalMesh(nx, 1., 3)

V = FunctionSpace(mesh, 'P', 1)

x = SpatialCoordinate(mesh)

# Define boundary condition
u_D = Expression("1./(pow(x[0]-y,2) + pow(eps, 2))", eps=1e-6, pi=pi, y=1., degree=2)
# plot(u_D)

def right(x, on_boundary):
    return near(x[0], 1.)
def left(x, on_boundary):
    return near(x[0], 1.)

bc = DirichletBC(V, Constant(0), right)

# Define initial value
u_n = interpolate(u_D, V)
C = assemble(u_n*dx)
u_n.assign(u_n/C)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

a_form = (u*v + (dl/theta)*(x[0]*x[0]-1)*u.dx(0)*v.dx(0))*dx
L_form = u_n*v*dx

# Time-stepping
u = Function(V)
l = 0
for n in range(num_steps):

    # Update current time
    l += dl

    # Compute solution
    solve(a_form == L_form, u, bc)
    
    vtkfile = File('foquesol.pvd')
    vtkfile << u

    # Plot solution
    plot(u)
    #plot(mesh)

    # Update previous solution
    u_n.assign(u)
    print(assemble(u*dx))
   

# Hold plot
plt.xlabel('$\\eta$')
plt.ylabel('$P(\\eta,L)$')
plt.legend(['$P(\\eta,L_0)$','$P(\\eta,L_1)$'], loc='upper right')
plt.show()

but when I plot this, I get plots for fixed L and varying \eta on the x axis and on the y axis I have P(L_{fixed},\eta). How can I plot P for fixed \eta and vary L on the x axis instead? Would I have to discretize the \eta variable instead of L? I hope not…Thanks in advance.

2

Hi,
you can just evaluate u(eta) at each time step and store it in a list/array.

1 Like