Ask Your Question
2

How do I perform a change of variables for a pde

asked 2011-11-03 11:02:07 -0500

rtrwalker gravatar image

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 flag offensive close merge delete

1 answer

Sort by » oldest newest most voted
4

answered 2011-11-14 13:59:45 -0500

rtrwalker gravatar image

updated 2011-11-14 14:02:39 -0500

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)
edit flag offensive delete link more

Comments

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

DSM gravatar imageDSM ( 2011-11-16 03:03:55 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2011-11-03 11:02:07 -0500

Seen: 1,020 times

Last updated: Nov 14 '11