<- previous index next ->
If you have code that works for double, extending to complex
may be easy.
Python evaluating a polynomial on "real" coefficients
test_polyfit.py polyfit in Python, needs numpy
from numpy import array
from numpy import polyfit
from numpy import polyval
import pylab
def f(x):
return 5.0 + 4.0*x + 3.0*x*x + 2.0*x*x*x + 1.0*x*x*x*x
print "test_polyfit.py a x,y,5) "
print "fit 5.0 + 4.0*x + 3.0*x*x + 2.0*x*x*x + 1.0*x*x*x*x"
xx = [0.0 for i in range(5)]
yy = [0.0 for i in range(5)]
for i in range(5):
xx[i] = 0.1*(i+1)
yy[i] = f(xx[i]) # function we will fit
print "i=",
print i,
print " ,xx=",
print xx[i],
print " ,yy=",
print yy[i]
print "p=polyfit(xx,yy,4) for 5 points"
p=polyfit(xx,yy,4) # 5 input values, 4th order
print "polyfit coefficients"
print p
print "backwards? largest power first, expected 5, 4, 3, 2, 1"
print " "
print "polyval values of fit"
yy1=polyval(p,xx)
print yy1
print "should be values above"
output
test_polyfit.py a x,y,5)
fit 5.0 + 4.0*x + 3.0*x*x + 2.0*x*x*x + 1.0*x*x*x*x
i= 0 ,xx= 0.1 ,yy= 5.4321
i= 1 ,xx= 0.2 ,yy= 5.9376
i= 2 ,xx= 0.3 ,yy= 6.5321
i= 3 ,xx= 0.4 ,yy= 7.2336
i= 4 ,xx= 0.5 ,yy= 8.0625
p=polyfit(xx,yy,4) for 5 points
polyfit coefficients
[ 1. 2. 3. 4. 5.]
backwards? largest power first, expected 5, 4, 3, 2, 1
polyval values of fit
[ 5.4321 5.9376 6.5321 7.2336 8.0625]
should be values above
Python evaluating a polynomial on "complex" coefficients
test_cxpolyfit.py complex polyfit in Python, needs numpy
from numpy import array
from numpy import polyfit
from numpy import polyval
import pylab
def f(x): # changed to some complex coefficients
return (5.0+1.0j) + (4.0+1.1j)*x + (3.0+1.2j)*x*x + 2.0*x*x*x + 1.0*x*x*x*x
print "test_cxpolyfit.py a x,y,5) using complex numbers "
print "fit (5.0+1.0j) + (4.0+1.1j)*x + (3.0+1.2j)*x*x + 2.0*x*x*x + 1.0*x*x*x*x"
xx = [0.0 for i in range(5)]
yy = [0.0 for i in range(5)]
for i in range(5):
xx[i] = 0.1*(i+1)
yy[i] = f(xx[i]) # now has complex coefficients
print "i=",
print i,
print " ,xx=",
print xx[i],
print " ,yy=",
print yy[i]
print "p=polyfit(xx,yy,4) for 5 points"
p=polyfit(xx,yy,4) # 5 input values, 4th order
print "polyfit coefficients"
print p
print "backwards? largest power first, expected about 5, 4, 3, 2, 1"
print " "
print "polyval values of fit"
yy1=polyval(p,xx)
print yy1
print "should be values above"
output
test_cxpolyfit.py a x,y,5) using complex numbers
fit (5.0+1.0j) + (4.0+1.1j)*x + (3.0+1.2j)*x*x + 2.0*x*x*x + 1.0*x*x*x*x
i= 0 ,xx= 0.1 ,yy= (5.4321+1.122j)
i= 1 ,xx= 0.2 ,yy= (5.9376+1.268j)
i= 2 ,xx= 0.3 ,yy= (6.5321+1.438j)
i= 3 ,xx= 0.4 ,yy= (7.2336+1.632j)
i= 4 ,xx= 0.5 ,yy= (8.0625+1.85j)
p=polyfit(xx,yy,4) for 5 points
polyfit coefficients
[ 1. +1.45680467e-12j 2. -1.77472663e-12j 3. +1.20000000e+00j
4. +1.10000000e+00j 5. +1.00000000e+00j]
backwards? largest power first, expected about 5, 4, 3, 2, 1
polyval values of fit
[ 5.4321+1.122j 5.9376+1.268j 6.5321+1.438j 7.2336+1.632j 8.0625+1.85j ]
should be values above
No difference in polyfit or polyval, they accept real or complex.
similar files to above
test_polyfit.py
test_polyfit_py.out
test_cxpolyfit.py
test_cxpolyfit_py.out
"C" language polyval on double
double polyval(int order, double p[], double x) // return y = p(x)
{ // using Horner's rule
double y;
int i;
y = p[order]*x; // p one larger than order
for(i=order-1; i>0; i--)
y = (p[i]+y)*x;
y = p[0]+y;
return y;
} // end polyval
"C" language polyval on complex
easy editing a*b becomes cxmul(a,b)
a+b becomes cxadd(a,b)
double becomes complex
// #include "complex.h" compile with complex.c
complex cxpolyval(int order, complex p[], complex x) // return y = p(x)
{ // using Horner's rule
complex y;
int i;
y = cxmul(p[order],x); // p one larger than order
for(i=order-1; i>0; i--)
y = cxmul(cxadd(p[i],y),x);
y = cxadd(p[0],y);
return y;
} // end cxpolyval
Java polyval on real
// polyval.java with simple test
class polyval
{
polyval(){}
double polyval(int order, double p[], double x) // return y = p(x)
{ // using Horner's rule
double y;
y = p[order]*x; // p one larger than order
for(int i=order-1; i>0; i--)
y = (p[i]+y)*x;
y = p[0]+y;
return y;
} // end polyval
public static void main (String[] args)
{
double y;
double x = 2.0;
double p[] = {1.0, 2.0, 3.0};
polyval c = new polyval();
y = c.polyval(2, p, x);
System.out.println("polyval.java on real");
System.out.println("p = {"+p[0]+","+p[1]+","+p[2]+"}");
System.out.println("y=p(x) x="+x+", y="+y);
} // end main
} // end polyval.java
polyval.java real output
polyval.java on real
p = {1.0,2.0,3.0}
y=p(x) x=2.0, y=17.0
Java polyval on complex
easy editing a*b becomes a.multiply(b)
a+b becomes a.add(b)
double becomes Complex
// cxpolyval.java with simple test on Complex
// needs class Complex.class that defines type and functions
class cxpolyval
{
cxpolyval(){}
Complex cxpolyval(int order, Complex p[], Complex x) // return y = p(x)
{ // using Horner's rule
Complex y;
y = p[order].multiply(x); // p one larger than order
for(int i=order-1; i>0; i--)
y = (p[i].add(y)).multiply(x);
y = p[0].add(y);
return y;
} // end cxpolyval
public static void main (String[] args)
{
Complex y;
Complex x = new Complex(2.0, 0.1);
Complex p[] = {new Complex(1.0,0.01), new Complex(2.0,0.02),
new Complex(3.0)};
cxpolyval c = new cxpolyval();
y = c.cxpolyval(2, p, x);
System.out.println("cxpolyval.java on Complex");
System.out.println("p = {"+p[0]+","+p[1]+","+p[2]+"}");
System.out.println("y=p(x) x="+x+", y="+y);
} // end main
} // end cxpolyval.java
cxpolyval.java complex output
cxpolyval.java on Complex
p = {(1.0,0.01),(2.0,0.02),(3.0,0.0)}
y=p(x) x=(2.0,0.1), y=(16.968,1.4500000000000002)
The notation (1.0,2.0) comes from Fortran where the
syntax was 1.0 real, 2.0 imaginary.
many Java files built for Complex
Complex.java
ComplexMatrix.java
ComplexPolynomial.java
Even least square fit can be performed on complex values
and the polynomial may have powers of all variables
* Y_actual X1 X2
* -------- ----- -----
* 32.5+i21.2 1.0+i1.0 2.5+i2.0
* 7.2+i5.3 2.0+i1.1 2.5
* 6.9+i1.1 3.0+i1.2 2.7
* 22.4+i0.4 2.2+i2.2 2.1+i2.1
* 10.4+i2.3 1.5+i3.4 2.0
* 11.3+i4.0 1.6+i1.0 2.0
* ... ... ... need at least as many as c's
*
* Find coefficients c0, c1, c2, ... such that
* Y_approximate = c0 + c1*X1 + c2*X2 + c3*X1^2 + c4*X1*X2 + c5*X2^2 +
* c6*X1^3 + c7*X1^2*X2 + c8*X1*X2^2 + c9*X2^3 + ...
* and such that the sum of (Y_actual - Y_approximate) squared is minimized.
* The same simultaneous equations are built, be sure to use a simultaneous
* equation solver where all arithmetic was changed to complex.
some examples with multiple variables to higher powers
(Sorry, modifying to complex is left as an exercise to the student.)
least_square_fit_2d.c
least_square_fit_3d.c
least_square_fit_4d.c
least_square_fit_4d.java
least_square_fit_4d_java.out
<- previous index next ->
-
CMSC 455 home page
-
Syllabus - class dates and subjects, homework dates, reading assignments
-
Homework assignments - the details
-
Projects -
-
Partial Lecture Notes, one per WEB page
-
Partial Lecture Notes, one big page for printing
-
Downloadable samples, source and executables
-
Some brief notes on Matlab
-
Some brief notes on Python
-
Some brief notes on Fortran 95
-
Some brief notes on Ada 95
-
An Ada math library (gnatmath95)
-
Finite difference approximations for derivatives
-
MATLAB examples, some ODE, some PDE
-
parallel threads examples
-
Reference pages on Taylor series, identities,
coordinate systems, differential operators
-
selected news related to numerical computation