Hi,
I'm trying to use sage to automatically derive the adjoint PDE
associated with the one-dimensional Burgers' equation. Suppose you
have an equation F(u) = 0, where u is the function you're trying to
solve for. Suppose you are interested in the value of some functional
J(u). The associated adjoint equation is formed by constructing the
Lagrangian:
L(u, lambda) = J(u) - <lambda, F(u)>
where <.> means the inner product.
The adjoint PDE is defined by diff(L, u) = 0. This is what I'm trying
to get sage to compute.
Unfortunately, I can't quite figure out the right construction:
# I am attempting to derive the adjoint PDE for the one-dimensional
Burgers
equation
# posed as the following:
# diff(u, t) + u * diff(u, x) - diff(u, x, 2)) = 0 on [a,
b] x [0, T]
# u(x, t=0) = u0 on [a,
b]
# u(a, t) = 0 on [0,
T]
# u(b, t) = 0 on [0,
T]
# To do this, I am attempting to form the Lagrangian
# L(u, rho) = J(u) - <rho, diff(u, t) + u * diff(u, x) - diff(u, x,
2))>_ST - <eta, u(x=a, t)>_T - <zeta, u(x=b), t)>_T - <chi, u(x,
t=0)>_S
# where _SP means inner product over spacetime, _S over space, and _T
over time.
# The adjoint equation is defined as
# diff(L, u) = 0
# which is what I would like to compute.
x = var('x') # space
t = var('t') # time
u = function('u', x, t) # velocity, to be solved for in the
Burgers equation
a = var('a') # the left-hand space limit
b = var('b') # the right-hand space limit
assume(a < b)
T = var('T') # the upper time limit
assume(T > 0)
def ip_st(u, v):
# inner product over spacetime
return integral(integral(u*v, x, a, b), t, 0, T)
def ip_t(u, v):
# inner product over time
return integral(u*v, t, 0, T)
def ip_s(u, v):
# inner product over space
return integral(u*v, x, a, b)
rho = function('rho', x, t) # the main adjoint variable
eta = function('eta', t) # the lagrange multiplier enforcing the
left boundary condition
zeta = function('zeta', t) # the lagrange multiplier enforcing the
right boundary condition
chi = function('chi', x) # the lagrange multiplier enforcing the
initial condition
# In the Lagrangian, we'll leave out the J(u), as I know what to do
with it
# Here is where it starts to go wrong.
L = -ip_st(rho, diff(u, t) + u * diff(u, x) - diff(u, x, 2)) -
ip_t(eta, u.subs(x=a)) - ip_t(zeta, u.subs(x=b)) - ip_s(chi,
u.subs(t=0))
# I am confused as to what to do now. L isn't a function of u, sage
says:
# sage: u in L.variables()
# False
# sage: print diff(L, u)
# ...
# TypeError: argument symb must be a symbol
eps = var('epsilon') # a small number
assume(eps > 0)
du = function('du', x, t) # a perturbation
# sage: print L.substitute_function(u, u + eps*du) - L
# 0
So I can't diff(L, u); it gives a TypeError. And when I try to
substitute_function(), it gives the same thing back again.
I tried many variants of this, making everything vars instead of
functions, but I couldn't hit upon the right trick.
Any ideas?
Thanks,
John
--
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
