Application Center - Maplesoft
  Submit your work
Application Center Applications Conversion linear partial differential operators with constant coefficients to canonical form Preview

App Preview:

Conversion linear partial differential operators with constant coefficients to canonical form

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application

 

pdetocan1.mws

Conversion linear partial differential operators
with constant coefficients to canonical form

by Aleksas Domarkas

Vilnius University, Faculty of Mathematics and Informatics,

Naugarduko 24, Vilnius, Lithuania

aleksas@ieva.mif.vu.lt

NOTE: In this session we find change of variables which reduct arbitrary linear second-order partial differential operator with constant coefficients to canonical form.  Program sqsum  convert quadratic form to sum squares.

Introduction

>    restart;

Please  input  number of examples  k (1..13) and   select    Edit->Execute->Worksheet

>    k:=11;

k := 11

>    with(linalg): 

Warning, the protected names norm and trace have been redefined and unprotected

Linear operators

>    L1:=D[1,1]+2*D[1,2]-2*D[1,3]+2*D[2,2]+6*D[3,3]:

>    L2:=D[1,1]+2*D[1,2]-2*D[1,3]+2*D[2,2]+2*D[3,3]:

>    L3:=D[1,1]+2*D[1,2]+2*D[2,2]+2*D[2,3]+2*D[2,4]+2*D[3,3]+3*D[4,4]:

>    L4:=D[1,2]-D[1,4]+D[3,3]-2*D[3,4]+2*D[4,4]:

>    L5:=D[1,1]+2*D[1,3]-2*D[1,4]+D[2,2]+2*D[2,3]+2*D[2,4]+2*D[3,3]+2*D[4,4]:

>    L6:=D[1,2]+D[1,3]+D[1,4]+D[3,4]:

>    L7:=D[1,1]+2*D[1,2]-2*D[1,3]-4*D[2,3]+2*D[2,4]+D[3,3]:

>    L8:=D[1,1]+2*D[1,2]-3*D[2,2]+2*D[1]+6*D[2]:

>    L9:=D[1,1]+4*D[1,2]+5*D[2,2]+D[1]+2*D[2]:

>    L10:=D[1,2]+D[1,3]+D[2,3]-D[1]+D[2]:

>    L11:=D[1,1]+4*D[1,2]+4*D[2,2]+6*D[1,3]+12*D[2,3]+9*D[3,3]-2*D[1]-4*D[2]-6*D[3]:

>    L12:=D[1,2]-2*D[1,3]+D[2,3]+D[1]+1/2*D[2]:

>   

>    L13:=D[1,2]:

>    L:=L||k;

L := D[1,1]+4*D[1,2]+4*D[2,2]+6*D[1,3]+12*D[2,3]+9*D[3,3]-2*D[1]-4*D[2]-6*D[3]

>   

Conversion quadratic form to sum squares

program sqsum

>    sqsum := proc(f::quadratic)
option `Copyright Aleksas Domarkas, 2000`;
local i, l, n, x, J, S, K, F, kk;
    if ldegree(f) <> 2 then error "f is not quadratic form"
    end if;
    S := f;
    K := 0;
    indets(f);
    x := convert(%, list);
    n := nops(x);
    while S <> 0 do
        while has(S, {seq(x[i]^2, i = 1 .. n)}) do for i to n
            do
                if has(S, x[i]^2) then
                    K := K +
                        1/4*diff(S, x[i])^2/coeff(S, x[i]^2);
                    S := expand(Q - K)
                end if
            end do
        end do;
        if S <> 0 then
            if type(S, `+`) then op(1, S) else S end if;
            indets(%);
            l := [coeff(coeff(%%, %[1]), %[2]), %[1], %[2]];
            K := K
                 + 1/4*(diff(S, l[2]) + diff(S, l[3]))^2/l[1]
                 - 1/4*(diff(S, l[2]) - diff(S, l[3]))^2/l[1]
                ;
            S := expand(f - K)
        end if
    end do;
    K := map(simplify, K);
    RETURN(K)
end proc;

>   

sqsum := proc (f::quadratic) local i, l, n, x, J, S, K, F, kk; option `Copyright Aleksas Domarkas, 2000`; if ldegree(f) <> 2 then error

examples using sqsum

Example 1

>    Q:=3*x1^2+3*x2^2+3*x3^2-2*x1*x2-2*x1*x3+2*x2*x3;

Q := 3*x1^2+3*x2^2+3*x3^2-2*x1*x2-2*x1*x3+2*x2*x3

>    sqsum(Q);simplify(%-Q);

1/3*(3*x1-x2-x3)^2+1/6*(4*x2+x3)^2+5/2*x3^2

0

Example 2

>    Q:=x*y+y*z+x*z;

Q := x*y+y*z+x*z

>    sqsum(Q);simplify(%-Q);

1/4*(y+2*z+x)^2-1/4*(-y+x)^2-z^2

0

Example 3

>    Q:=sum(sum(x[i]*x[j],j=1..i-1),i=1..8);

Q := x[2]*x[1]+x[3]*x[1]+x[3]*x[2]+x[4]*x[1]+x[4]*x[2]+x[4]*x[3]+x[5]*x[1]+x[5]*x[2]+x[5]*x[3]+x[5]*x[4]+x[6]*x[1]+x[6]*x[2]+x[6]*x[3]+x[6]*x[4]+x[6]*x[5]+x[7]*x[1]+x[7]*x[2]+x[7]*x[3]+x[7]*x[4]+x[7]*x...
Q := x[2]*x[1]+x[3]*x[1]+x[3]*x[2]+x[4]*x[1]+x[4]*x[2]+x[4]*x[3]+x[5]*x[1]+x[5]*x[2]+x[5]*x[3]+x[5]*x[4]+x[6]*x[1]+x[6]*x[2]+x[6]*x[3]+x[6]*x[4]+x[6]*x[5]+x[7]*x[1]+x[7]*x[2]+x[7]*x[3]+x[7]*x[4]+x[7]*x...
Q := x[2]*x[1]+x[3]*x[1]+x[3]*x[2]+x[4]*x[1]+x[4]*x[2]+x[4]*x[3]+x[5]*x[1]+x[5]*x[2]+x[5]*x[3]+x[5]*x[4]+x[6]*x[1]+x[6]*x[2]+x[6]*x[3]+x[6]*x[4]+x[6]*x[5]+x[7]*x[1]+x[7]*x[2]+x[7]*x[3]+x[7]*x[4]+x[7]*x...

>    sqsum(Q);simplify(%-Q);

1/4*(x[2]+2*x[3]+2*x[4]+2*x[5]+2*x[6]+2*x[7]+2*x[8]+x[1])^2-1/4*(-x[2]+x[1])^2-1/4*(x[3]+x[5]+x[6]+x[7]+x[8]+2*x[4])^2-1/12*(x[3]+x[6]+x[7]+x[8]+3*x[5])^2-1/24*(x[6]+x[7]+x[8]+4*x[3])^2-1/40*(x[7]+x[8]...
1/4*(x[2]+2*x[3]+2*x[4]+2*x[5]+2*x[6]+2*x[7]+2*x[8]+x[1])^2-1/4*(-x[2]+x[1])^2-1/4*(x[3]+x[5]+x[6]+x[7]+x[8]+2*x[4])^2-1/12*(x[3]+x[6]+x[7]+x[8]+3*x[5])^2-1/24*(x[6]+x[7]+x[8]+4*x[3])^2-1/40*(x[7]+x[8]...
1/4*(x[2]+2*x[3]+2*x[4]+2*x[5]+2*x[6]+2*x[7]+2*x[8]+x[1])^2-1/4*(-x[2]+x[1])^2-1/4*(x[3]+x[5]+x[6]+x[7]+x[8]+2*x[4])^2-1/12*(x[3]+x[6]+x[7]+x[8]+3*x[5])^2-1/24*(x[6]+x[7]+x[8]+4*x[3])^2-1/40*(x[7]+x[8]...

0

Example 4

>    n:=5:A:=randmatrix(n,n,symmetric):X:=vector(n,[x||(1..n)]):

>    Q:=expand(evalm(transpose(X)&*A&*X));

Q := -110*x1*x2-70*x1*x3+194*x2*x3-85*x1^2-37*x2^2+50*x3^2+158*x1*x4+114*x1*x5+112*x2*x4-118*x2*x5+98*x3*x4+90*x3*x5+63*x4^2-16*x4*x5-93*x5^2
Q := -110*x1*x2-70*x1*x3+194*x2*x3-85*x1^2-37*x2^2+50*x3^2+158*x1*x4+114*x1*x5+112*x2*x4-118*x2*x5+98*x3*x4+90*x3*x5+63*x4^2-16*x4*x5-93*x5^2

>    sqsum(Q);simplify(%-Q);

-1/85*(55*x2+35*x3+85*x1-79*x4-57*x5)^2-1/408*(-2034*x3+24*x2-83*x4+1630*x5)^2+1/163272*(40818*x3+1721*x4-32418*x5)^2+1/22520430646740*(55172793*x4+22487254*x5)^2-311683396/165518379*x5^2
-1/85*(55*x2+35*x3+85*x1-79*x4-57*x5)^2-1/408*(-2034*x3+24*x2-83*x4+1630*x5)^2+1/163272*(40818*x3+1721*x4-32418*x5)^2+1/22520430646740*(55172793*x4+22487254*x5)^2-311683396/165518379*x5^2
-1/85*(55*x2+35*x3+85*x1-79*x4-57*x5)^2-1/408*(-2034*x3+24*x2-83*x4+1630*x5)^2+1/163272*(40818*x3+1721*x4-32418*x5)^2+1/22520430646740*(55172793*x4+22487254*x5)^2-311683396/165518379*x5^2

0

Example 5

>    sqsum(x*y+1); 

Error, (in sqsum) f is not quadratic form

>    sqsum(x*y);

1/4*(y+x)^2-1/4*(-y+x)^2

>   

Conversion differential operator  to canonical form

>    L;

D[1,1]+4*D[1,2]+4*D[2,2]+6*D[1,3]+12*D[2,3]+9*D[3,3]-2*D[1]-4*D[2]-6*D[3]

>    indets(L,indexed);

{D[1,1], D[1,2], D[1,3], D[2,2], D[3,3], D[2,3], D[1], D[2], D[3]}

>    map(op,%);

{1, 2, 3}

>    n:=max(op(%));

n := 3

>    Q:=subs(seq(seq(D[i,j]=x[i]*x[j],i=1..n),j=1..n),seq(D[j]=0,j=1..n),L);

Q := x[1]^2+4*x[2]*x[1]+4*x[2]^2+6*x[3]*x[1]+12*x[3]*x[2]+9*x[3]^2

>    A:=hessian(Q/2,[seq(x[j],j=1..n)]);

A := matrix([[1, 2, 3], [2, 4, 6], [3, 6, 9]])

>    C:=[seq(coeff(L,D[j]),j=1..n)];

C := [-2, -4, -6]

>    X :=[seq(x[i],i=1..n)];

X := [x[1], x[2], x[3]]

>    K:=sqsum(Q);

K := (x[1]+2*x[2]+3*x[3])^2

>    F:=x->subs(I=1,sqrt(x,symbolic)):

>    if type(K,`+`) then map(F,[op(K)]) else [F(K)] end if;

[x[1]+2*x[2]+3*x[3]]

>    kk := [op(%),seq(X[i],i=nops(%)+1..n)];

kk := [x[1]+2*x[2]+3*x[3], x[2], x[3]]

>    jacobian(kk,X);

matrix([[1, 2, 3], [0, 1, 0], [0, 0, 1]])

>    B :=transpose(evalm(( jacobian(kk,X))^(-1)));

B := matrix([[1, 0, 0], [-2, 1, 0], [-3, 0, 1]])

>    itr :=factor({seq(y[i]=multiply(B,X)[i],i=1..n)});

itr := {y[1] = x[1], y[3] = -3*x[1]+x[3], y[2] = -2*x[1]+x[2]}

>    An :=evalm( B&*A&*transpose(B));

An := matrix([[1, 0, 0], [0, 0, 0], [0, 0, 0]])

>    Bn:=[seq(sum(B[k,i]*C[i],i=1..n),k=1..n)];

Bn := [-2, 0, 0]

>    L_can :=sum(sum(An[i,j]*D[i,j],i=1..n),j=1..n)+sum(Bn[i]*D[i],i=1..n);

L_can := D[1,1]-2*D[1]

Example

>    cat(k,` EXAMPLE `);'L'=L; 'L_can'=L_can;'itr'=itr;

`11 EXAMPLE `

L = D[1,1]+4*D[1,2]+4*D[2,2]+6*D[1,3]+12*D[2,3]+9*D[3,3]-2*D[1]-4*D[2]-6*D[3]

L_can = D[1,1]-2*D[1]

itr = {y[1] = x[1], y[3] = -3*x[1]+x[3], y[2] = -2*x[1]+x[2]}

>   

Conversion to diff and checking    with PDEtools[dchange]

>    L:=convert(L(u)(seq(x[k],k=1..n)),diff);

L := diff(u(x[1],x[2],x[3]),`$`(x[1],2))+4*diff(u(x[1],x[2],x[3]),x[1],x[2])+4*diff(u(x[1],x[2],x[3]),`$`(x[2],2))+6*diff(u(x[1],x[2],x[3]),x[1],x[3])+12*diff(u(x[1],x[2],x[3]),x[2],x[3])+9*diff(u(x[1]...
L := diff(u(x[1],x[2],x[3]),`$`(x[1],2))+4*diff(u(x[1],x[2],x[3]),x[1],x[2])+4*diff(u(x[1],x[2],x[3]),`$`(x[2],2))+6*diff(u(x[1],x[2],x[3]),x[1],x[3])+12*diff(u(x[1],x[2],x[3]),x[2],x[3])+9*diff(u(x[1]...
L := diff(u(x[1],x[2],x[3]),`$`(x[1],2))+4*diff(u(x[1],x[2],x[3]),x[1],x[2])+4*diff(u(x[1],x[2],x[3]),`$`(x[2],2))+6*diff(u(x[1],x[2],x[3]),x[1],x[3])+12*diff(u(x[1],x[2],x[3]),x[2],x[3])+9*diff(u(x[1]...

>    convert(L_can(u)(seq(y[k],k=1..n)),diff);

diff(u(y[1],y[2],y[3]),`$`(y[1],2))-2*diff(u(y[1],y[2],y[3]),y[1])

Checking  with PDEtools[dchange]:

>    tr:=solve(itr,{seq(x[k],k=1..n)}):

>    PDEtools[dchange](tr,L,itr,[seq(y[k],k=1..n)],simplify);

diff(u(y[1],y[2],y[3]),`$`(y[1],2))-2*diff(u(y[1],y[2],y[3]),y[1])

>   

While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.

Back to contents