Project 1 will be to implement a finite element method for
a two-point boundary-value problem. It will

have several parts.

**Warmup: Solving quadratic equations**

The quadratic formula says that the solutions of ax ^{2} + bx + c = 0 are
given by

If b^{2} − 4ac < 0 then √
is imaginary, so there is no problem with round-off error.

If b^{2} − 4ac > 0 then cancellation can occur in
if b > 0 and in
if b < 0.

Thus, if b > 0 one would want to use the computations

Similarly if b < 0.

Write a function (quadratic-solver a b c) that returns a list of the two roots
of ax ^{2} + bx + c = 0 as

accurately as possible. Test your code on the following problems:

euler-126% gsi

Gambit v4.1.2

> (sqrt -1)

+i

> (sqrt +i)

.7071067811865476+.7071067811865475i

> (load "quadratic-solver")

"/export/users/lucier/programs/615project/2007/project-1/quadratic-solver.scm"

> (quadratic-solver 1 2 5)

(-1+2i -1-2i)

> (quadratic-solver 1 -1 1)

(1/2+.8660254037844386i 1/2-.8660254037844386i)

> (quadratic-solver 1 2 -1)

(-2.414213562373095 .4142135623730951)

> (quadratic-solver 4 1 1)

(-1/8+.4841229182759271i -1/8-.4841229182759271i)

> (quadratic-solver 4 4 1)

(-1/2 -1/2)

> (quadratic-solver 4 0 1)

(+1/2i -1/2i)

> (quadratic-solver 0 0 1)

*** ERROR IN (console)@11.1 -- not a quadratic: 0 0 1

1>

> (quadratic-solver 4 0 -1)

(-1/2 1/2)

> (quadratic-solver 1 3138428376721 1)

(-3.138428376721e12 -3.186308177103568e-13)

>

*** EOF again to exit

**Meroon**

Standard Scheme (so-called R5RS Scheme, which Gambit implements) does not have
an object system.

We use an object system provided by the software package Meroon.

To use Gambit, you need to have /pkgs/Gambit-C/current/bin/ in your path. The
Gambit interpreter

is called gsi and the Gambit compiler is called gsc.

To have Gambit load Meroon automatically, just call gsi++ or gsc++.

Our system has two differences with standard Meroon:

(1) In standard Meroon, keywords begin with a colon; in our Meroon keywords end
with a colon:

(define-class Polynomial Object

((= variable immutable:)

(= terms immutable:)))

(2) In standard Meroon, so-called setters begin with set- and end with !. In our
Meroon, setters end

with -set!:

euler-130% gsi++

[ Meroon V3 Paques2001+1 $Revision: 1.1 $ ]

Gambit v4.1.2

> (define-class Point Object (x y))

Point

> (define p (make-Point 0 1))

> (unveil p)

(a Point <------------- [Id: 1]

x: 0

y: 1 end Point)

#t

> (Point-x-set! p 1)

#<meroon #2>

> (unveil p)

(a Point <------------- [Id: 1]

x: 1

y: 1 end Point)

#t

>

** Numerical Integration **

This first part will be about numerical integration (quadrature rules).

The Gauss-Lobatto quadrature rules with n points have the form

Here are the zeros of
the degree n − 2 orthogonal polynomial over [−1, 1] with the weight

If we define

then has degree n − 1
and satisfies

The weights satisfy

So, the first part of the project is to write code to
manipulate polynomials.

and modify it to use Meroon’s framework of classes/objects and generics/methods.

We’ll define a polynomial class:

(define-class Polynomial Object

((= variable immutable:)

(= terms immutable:)))

and a way to check whether two Polynomial variables are the same:

(define (Polynomial-variable= var1 var2)

(eq? var1 var2))

The terms of a polynomial is just a list of nonzero terms, in decreasing order
by degree (unfortunately called

“order” at that web page), so we need some code to manipulate terms and lists of
terms:

;;; a term is a pair (coeff order) (order should really be degree , but ...)

;;; (Polynomial-terms p) is a list of terms in decreasing orders.

;; ope ration on terms and term-lists

(define (adjoin-term term term-list)

(if (=zero? (term-coeff term))

term-list

(cons term term-list)))

(define (the-empty-termlist)

’())

(define (first-term term-list)

(car term-list))

(define (rest-terms term-list)

(cdr term-list))

(define (empty-termlist? term-list)

(null? term-list))

(define (make-term order coeff)

(list order coeff))

(define (term-order term)

(car term))

(define (term-coeff term)

(cadr term))

The web page has code for adding two polynomials . Putting it into our terms we
define a generic function

add that should work for everything, and we start with it working with numbers:

(define-generic (add (x) y)

(if (and (number? x)

(number? y))

(+ x y)

(error "add: This generic is not defined on these objects: " x y)))

and then we define a method for adding Polynomials:

(define-method (add (p_1 Polynomial) p_2)

(cond ((number? p_2)

(add p1 (number->Polynomial p2 (Polynomial-variable p1))))

((and (Polynomial? p_2)

(Polynomial-variable= (Polynomial-variable p_1)

(Polynomial-variable p_2)))

(instantiate Polynomial

variable: (Polynomial-variable p_1)

terms: (add-terms (Polynomial-terms p_1)

(Polynomial-terms p_2))))

(else

(error "add: p_2 is neither a number nor a polynomial with the same variable as

p_1 " p_1 p_2))))

This method is called only when p_1 is a polynomial; if p_2 is a number, it
converts p _2 to a Polynomial

with the same variable as p_1 and calls add again with both arguments now a
Polynomial.

The web page has code for add-terms:

(define (add-terms l1 l2)

(cond ((empty-termlist? l1) l2)

((empty-termlist? l2) l1)

(else

(let ((t1 (first-term l1))

(t2 (first-term l2)))

(cond ((> (term-order t1)

(term-order t2))

(adjoin-term t1

(add-terms (rest-terms l1) l2)))

((< (term-order t1)

(term-order t2))

(adjoin-term t1

(add-terms l1 (rest-terms l2))))

(else

(adjoin-term

(make-term (term-order t1)

(add (term-coeff t1)

(term-coeff t2)))

(add-terms (rest-terms l1)

(rest-terms l2)))))))))

So you need to define number->Polynomial, which takes two arguments.

You need to define a multiply generic that works with numbers by default, and a
method for multiply

that works on Polynomials; follow ths same pattern as for add. The web page has
the guts of the code:

(define (multiply-terms l1 l2)

(if (empty-termlist? l1)

(the-empty-termlist)

(add-terms (multiply-term-by-all-terms (first-term l1) l2)

(multiply-terms (rest-terms l1) l2))))

(define (multiply-term-by-all-terms t1 L)

(if (empty-termlist? L)

(the-empty-termlist)

(let ((t2 (first-term L)))

(adjoin-term

(make-term (+ (term-order t1)

(term-order t2))

(multiply (term-coeff t1)

(term-coeff t2)))

(multiply-term-by-all-terms t1 (rest-terms L))))))

So that’s pretty much the code that comes on the web page. Meroon defines a
generic function show

that we can specialize for Polynomials as such:

(define-method (show (p Polynomial) . stream)

(let ((port (if (null? stream)

(current-output-port)

(car stream))))

(if (=zero? p)

(display 0)

(show-terms (Polynomial-variable p)

(Polynomial-terms p)

port))

(newline port)))

(define (show-terms variable terms port)

(show-first-term variable (first-term terms) port)

(for-each (lambda (term)

(show-term variable term port))

(rest-terms terms)))

(define (show-first-term variable term port)

(let ((coeff (term-coeff term))

(order (term-order term)))

(display (list (if (and (= coeff 1)

(positive? order))

’()

coeff)

(cond ((zero? order)’())

((= order 1) variable)

(else

(list variable "^" order)))))))

(define (show-term variable term port)

(let ((coeff (term-coeff term))

(order (term-order term)))

(display (list (if (negative? coeff)

coeff

(list "+" coeff))

(cond ((zero? order)’())

((= order 1) variable)

(else

(list variable "^" order))))

port)))

It will probably help your debugging.

So, here are some problems.

(1) The above code uses a function =zero?. Define a generic function =zero? that
handles numbers.

Define a method that works with Polynomials.

(2) Define a generic function (negate (x)) that handles numbers by default.
Define a method for

negate that works with Polynomials. Use the generic negate to define a regular
function ( subtract

x y).

(3) Define a function (exponentiate x n) that uses multiply to exponentiate
anything that multiply

can multiply. Use the discussion of exponentiation on page

as your model.

(4) Define a function (variable->Polynomial x) that takes a symbol x and returns
a Polynomial that

represents the polynomial x, i.e., a single term with coefficient 1 and order 1.

(5) Define a generic function (evaluate f x) that evaluates the function f at x.
If f is a number,

assume that it means a function that constantly returns f. Define a method for
Polynomials.

If you’ve done the exercises until now, something like the following should
work.

;;; evaluation

(define-generic (evaluate (f) x)

(if (number? f)

f

(error "evaluate: unknown argument types " f x)))

(define-method (evaluate (p Polynomial) x)

(evaluate-terms (Polynomial-terms p) x))

(define (evaluate-terms terms x)

(if (empty-termlist? terms)

0

(add (evaluate-term (first-term terms) x)

(evaluate-terms (rest-terms terms) x))))

(define (evaluate-term term x)

(multiply (exponentiate x (term-order term))

(term-coeff term)))

Can you write a method that uses Horner’s rule for evaluating Polynomials in our
representation?