>> From: Richard Patterson "rpatters"
----------
From: Richard Patterson <rpatters>
Date: Mon, 17 Feb 2003 15:22:29 -0500 (EST)
To: "helmut.kahovec" "fwchapma"
"devore" "sylvester7"
Cc: "rpatters"
Subject: polynomials to operators
Dear Fred, Helmut, Chris, Francis and now Carl,
Below is what I came up with after your many helpful suggestions.
My particular application uses homogeneous polynomials, so is
easier for that reason.
I admit to having a bias toward code that I completely understand.
> restart;
> f:=(x,y)->a*x^2+b*x*y+c*y^2;
2 2
f := (x, y) -> a x + b x y + c y
> g:=(x,y)->p*x^3+q*x^2*y+r*x*y^2+s*y^3;
3 2 2 3
g := (x, y) -> p x + q x y + r x y + s y
Procedure to change the coefficients of polynomials to constant functions.
> Change:=proc(f) local C:
> for C in coeffs(f(x,y),[x,y]) do assign(C,unapply(C,f)) end do;
> end proc;
Change := proc(f)
local C;
for C in coeffs(f(x, y), [x, y]) do
assign(C, unapply(C, f))
end do
end proc
This leaves me uneasy as it affects f even though the procedure returns
nothing. However the Changed polynomials can still be multiplied,
evaluated, differentiated just like ordinary polynomials.
> Test. This is how I know Change does something to f. There is no
other reason for this evaluation.
> f(x,y)(3);
2 2
a(3) x(3) + b(3) x(3) y(3) + c(3) y(3)
> Change(f);
> f(x,y)(3);
2 2
a x(3) + b x(3) y(3) + c y(3)
> Procedure to create an operator from a homogeneous polynomial:
> Operator:=proc(f) local i,j,n,DD;
> n:=degree(f(x,y),{x,y});
> DD:=0;
> for i from 0 to n do
> DD:=DD+coeff(coeff(f(x,y),x,i),y,n-i)*D[seq(1,j=1..i),seq(2,j=1..n-i)]
> od;
> RETURN(DD);
> end;
Operator := proc(f)
local i, j, n, DD;
n := degree(f(x, y), {x, y});
DD := 0;
for i from 0 to n do DD := DD +
coeff(coeff(f(x, y), x, i), y, n - i)*
D[seq(1, j = 1 .. i), seq(2, j = 1 .. n - i)]
end do;
RETURN(DD)
end proc
> Test:
> Operator(f);
c D[2, 2] + b D[1, 2] + a D[1, 1]
> Operator(f)(f)(x,y);
2 2 2
2 c + b + 2 a
> Operator(f)(g);
c ((x, y) -> 2 r x + 6 s y) + b ((x, y) -> 2 q x + 2 r y)
+ a ((x, y) -> 6 p x + 2 q y)
> Operator(f)(g)(x,y);
c (2 r x + 6 s y) + b (2 q x + 2 r y) + a (6 p x + 2 q y)
> It is also possible to define the operator before using Change.
> Operator(g);
s D[2, 2, 2] + r D[1, 2, 2] + q D[1, 1, 2] + p D[1, 1, 1]
> Operator(g)(g)(x,y);
6 s(g)(x, y) s + 2 r(g)(x, y) r + 2 q(g)(x, y) q + 6 p(g)(x, y) p
> Change(g);
> Operator(g)(g)(x,y);
2 2 2 2
6 s + 2 r + 2 q + 6 p
> Operator(g)(f)(x,y);
0
This was expected, since degree of g is bigger than degree of f.
|
Engineering