# How do I perform a change of variables for a pde

How can I transform, step by step, the partial differential equation (pde) u.diff(z,z)k/gw == mvu.diff(t) into u.diff(Z,Z)=u.diff(Tv) where Z=z/H, Tv = cvt/H^2 and cv = k/(gwmv)? As you can see from my attempt below I can isolate cv; after that I have no idea how to do the change of variables. I can do this example by hand: let Z=z/H --> du/dz = du/dZ * dZ/dz --> du/dz = du/dZ *(1/H) ...etc. But Other examples may not be so easy (and I don't want to do it by hand!).

var('u,z,t,Z,H,Tv,k,gw,mv,cv')
u = function('u',z,t)
eq1 = u.diff(z,z)*k/gw == mv * u.diff(t)
print(eq1)
mycoeff=eq1.left_hand_side().coeff(u.diff(z,z),1)
print(mycoeff)
eq2 = eq1.divide_both_sides(mycoeff)
print(eq2)
mycoeff = eq2.right_hand_side().coeff(u.diff(t),1)
eq3 = cv == 1/mycoeff
print(eq3)
eq4 = solve(eq3,eq3.right_hand_side().default_variable(), solution_dict=true)
print(eq4)
eq5 = eq2.subs(eq4[0])
print(eq5)


which gives me the following output:

k*D[0, 0](u)(z, t)/gw == mv*D[1](u)(z, t)
k/gw
D[0, 0](u)(z, t) == gw*mv*D[1](u)(z, t)/k
cv == k/(gw*mv)
[{gw: k/(cv*mv)}]
D[0, 0](u)(z, t) == D[1](u)(z, t)/cv

edit retag close merge delete

Sort by » oldest newest most voted

With some insights from the following posts:

substitute expression instead of formal function symbol

How do I get an ordered list of a symbolic functions arguments?

Call a function with argument list in python (non-sage link)

I've come up with a function, transform_func_x(expr,f,rels), to change the dependent variables of an expression. Below is the function itself and below that how I've used it (all copy and pasted from a sage notebook). Any suggested improvements are welcome.

def transform_func_x(expr,f,rels):
r"""
Transform the dependent variables of a symbolic function

INPUT:

- expr - symbolic function.  Usually a relational expression
with differential operators such as a partial differential equation.

- f -  expression defining a function (e.g. previously defined
as f = function ('u', z, t))

- rels - list of relational expressions giving variable
transformations.  For example to change y(z, t) to y(Z, T_v), rels
could be [Z == z/H, T_v = c_v * t / H^2].  To change y(z, t) to y(Z, t),
rels could be [Z == z / H]

.. note::

1. Transformations will not be performed for any relationship where the
LHS is already one of the variables in f. For, example if
f = u(Z, t) and rels = [Z == z / H, T_v == c_v * t / H^2] then the
"Z == z / H" transformation will be ignored.
2. After applying this function to expr, you will probably want to
update f by applying transform_func_x to f itself
i.e transform_func_x(f,f,rels)
"""
#Limit#all#lines#to#a#maximum#of#79#characters#################################
y = f.operator()
x = f.operands()
nx = f.number_of_operands()

new_arg_list = x[:]
"""used to transform variables within f to functions
of those variables. e.g. u(z, t) --> u(z / H, t)"""
subs_dict1 = {}
"""used to change name of transformed variables within f.
e.g. u(z / H, t) --> u(Z, t)"""
subs_dict2 = {}
"""used to transform variables outside of f.
e.g. z *u(Z, t)--> z * H * u(Z, t)"""
for v in rels:
if v.left() not in x:
"""i.e. exclude transformations where LHS of
rels is already a variable of f"""
for i in range(nx):
# find which variable the transformation applies to
if x[i] in v.right().variables():
xi=i
break
new_arg_list[xi] = v.right()
subs_dict1[v.right()]=v.left()
subs_dict2[x[xi]]=solve(v,x[xi],solution_dict=true)[_sage_const_0 ][x[xi]]

# make dummy function for transforming variables of f
dummy_arg_list = [var("x{0:d}".format(i)) for i in range(nx)]
dummy_y(*dummy_arg_list) = function(str(y),*new_arg_list))

return (expr.substitute_function(y,dummy_y).subs(subs_dict1).
subs(subs_dict2))

reset
var('z,t,H,T_v,Z,c_v,k_v,m_v, g_w')
f=function('u',z,t)
pde = k_v*f.diff(z,z)/g_w == m_v*f.diff(t)

print('Variable transformations')

print('PDE with u(z,t):')
rels=[]
eq1=transform_func_x(pde,f,rels)
f=transform_func_x(f,f,rels)
print(eq1)
print('Term to isolate:')
myterm = f.diff(f.operands()[0],2)
print (myterm)
print('Coefficient of term on LHS:')
mycoeff=eq1.left_hand_side().coeff(myterm)
print(mycoeff ...
more

BTW, many thanks for coming back and posting your solution! Often people forget about a question and leave it hanging, but this helps everybody.

( 2011-11-16 03:03:55 -0600 )edit