Hi Kent,
On Nov 23, 2007 9:04 AM, kent-and <[EMAIL PROTECTED]> wrote:
>
>
> I have made a small script that computes the basis of a Lagragian
> finite element in 2D
> and would like to integrate its basis functions.
>
> Basically, its like this (this works)
>
> f = 1 - 3*x - 3*y + 2*x**2 + 2*y**2 + 4*x*y
> intf = integrate(f, (y, 0, 1))
> intff = integrate(intf, (x, 0, 1))
>
> When I compute the Lagrangian basis, the first basis function looks
> like f, when I print it out. However, there is obviously some
> difference between
> f and the first basis function, because when I try to integrate I
> get:
>
> Traceback (most recent call last):
> File "sympy_test.py", line 90, in <module>
> intff = integrate(intf, (x, 0, 1))
> File "/usr/lib/python2.5/site-packages/sympy/integrals/
> integrals.py", line 139, in integrate
> return integral.doit()
> File "/usr/lib/python2.5/site-packages/sympy/integrals/
> integrals.py", line 90, in doit
> antideriv = self._eval_integral(function, x)
> File "/usr/lib/python2.5/site-packages/sympy/integrals/
> integrals.py", line 120, in _eval_integral
> return risch_norman(f, x)
> File "/usr/lib/python2.5/site-packages/sympy/integrals/risch.py",
> line 182, in risch_norman
> if not f.has(x):
> File "/usr/lib/python2.5/site-packages/sympy/core/basic.py", line
> 387, in has
> p = Basic.sympify(patterns[0])
> File "/usr/lib/python2.5/site-packages/sympy/core/basic.py", line
> 251, in sympify
> raise NotImplementedError('matrix support')
> NotImplementedError: matrix support
>
>
>
>
>
> The complete script is as follows:
>
>
> from sympy import *
>
> x, y = symbols('xy')
>
> def integral_over_reference_square(f):
> intf = integrate(f, (y, 0, 1))
> intff = integrate(intf, (x, 0, 1))
> return intff
>
>
>
> def pol_space(order):
> sum = 0
> basis = []
> coeff = []
> for o1 in range(0,order+1):
> for o2 in range(0,order+1):
> if o1 + o2 <= order:
> aij = Symbol("a_%d_%d" % (o1,o2))
> sum += aij*pow(x, o1)*pow(y, o2)
> basis.append(pow(x, o1)*pow(y, o2))
> coeff.append(aij)
>
> return sum, coeff, basis
>
>
> def create_point_set(order):
> h = 1.0/(order)
> set = []
> for i in range(0, order+1):
> x = i*h
> for j in range(0, order+1):
> y = j*h
> if x + y <= 1.0 + h/2.0:
> set.append((x,y))
> return set
>
> def create_matrix(equations, coeffs):
> A = zeronm(len(equations), len(equations))
> i = 0; j = 0
> for j in range(0, len(coeffs)):
> c = coeffs[j]
> for i in range(0, len(equations)):
> e = equations[i]
> d, r = div(e, c)
> A[i,j] = d
>
> return A
>
>
>
>
> order = 2
> pol, coeffs, basis = pol_space(order)
> print "polynom ", pol
> print "polynom ", coeffs
> print "polynom ", basis
>
> points = create_point_set(order)
> print "points ", points
>
> equations = []
> for p in points:
> ex = pol.subs(x, p[0]).subs(y, p[1])
> equations.append(ex )
>
>
> A = create_matrix(equations, coeffs)
> Ainv = A.inv()
>
> b = eye(len(equations))
>
> x = Ainv*b
>
> print x
> basis = []
> for i in range(0,len(equations)):
> N = pol
> for j in range(0,len(coeffs)):
> N = N.subs(coeffs[j], x[j,i])
> basis.append(N)
>
>
>
> f = basis[0]
> print f
> intf = integrate(f, (y, 0, 1))
> print intf
> intff = integrate(intf, (x, 0, 1))
> print intff
That's a very nice application. And that's an interesting problem - this works:
In [1]: f = 0.16666 - x + 2*x**2
In [2]: integrate(f, (x, 0, 1))
Out[2]: 0.3333266666666666692512658678
But your script, which does the same, doesn't. I am investigating it
at the moment.
Ondrej
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sympy?hl=en
-~----------~----~----~----~------~----~------~--~---
