Application Center - Maplesoft

# Conversion linear partial differential operators with constant coefficients to canonical form

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

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;

 > 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;

 >

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;

 >

examples using sqsum

Example 1

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

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

Example 2

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

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

Example 3

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

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

Example 4

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

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

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

Example 5

 > sqsum(x*y+1);

```Error, (in sqsum) f is not quadratic form
```

 > sqsum(x*y);

 >

Conversion differential operator  to canonical form

 > L;

 > indets(L,indexed);

 > map(op,%);

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

 > 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);

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

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

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

 > K:=sqsum(Q);

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

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

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

 > jacobian(kk,X);

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

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

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

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

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

Example

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

 >

Conversion to diff and checking    with PDEtools[dchange]

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

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

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);

 >

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.