← Back to branch summary

~daniele-bigoni/tensortoolbox/tt-docs

« back to all changes in this revision

Viewing changes to Development/Examples/QuadratureRules/MaxVolPointSelection.py

  • Committer: Daniele Bigoni
  • Date: 2015-02-08 10:36:29 UTC
  • Revision ID: dabi@dtu.dk-20150208103629-61jivgkcgk4ksqr2
working on quadrature example

Show diffs side-by-side

added added

removed removed

Lines of Context:
Development/Examples/QuadratureRules/MaxVolPointSelection.py
30
30
#  Try to throw l points on a n-points Gauss rule
31
31
#  
32
32
 
33
 
TESTS = [1]
 
33
TESTS = [4]
34
34
 
35
35
 
36
36
if 1 in TESTS:
315
315
    plt.imshow( np.log10(np.abs(I_tmp)), interpolation='none')
316
316
    plt.colorbar()
317
317
    plt.show(block=False)
 
318
 
 
319
if 4 in TESTS:
 
320
    ####################################################
 
321
    # TEST 4
 
322
 
 
323
    xspan = [-1.,1.]
 
324
    l = 1
 
325
    N = 10000
 
326
    X = np.linspace(-1.,1.,N)
 
327
    tot_it = 5
 
328
 
 
329
    P = S1D.Poly1D(S1D.JACOBI,[0.,0.])
 
330
 
 
331
    (xG,wG) = P.GaussQuadrature(l,normed=True)
 
332
 
 
333
    x_list = [ xG ]
 
334
    for i in range(tot_it):
 
335
        # Compute moments of orthonormal polys
 
336
        npoint = (l+1)**(i+2) - 1
 
337
        (x,w) = P.GaussQuadrature(npoint,normed=True)
 
338
        v = P.GradVandermonde1D(x,npoint,0,norm=True)
 
339
        b = np.dot(v.T,w)
 
340
 
 
341
        # Throw additional points
 
342
        X  = np.hstack((x_list[-1],np.linspace(xspan[0],xspan[1],N)))
 
343
 
 
344
        # Construct Vandermonde matrix
 
345
        V_cp = P.GradVandermonde1D(X,npoint,0,norm=True)
 
346
 
 
347
        # Bias Gauss Points
 
348
        V_cp[:len(x_list[-1]),:] *= 10.
 
349
 
 
350
        # Run Maxvol
 
351
        (I, VsqInv, it) = DT.maxvol(V_cp, delta=1e-10,maxit=1000)
 
352
 
 
353
        # Compute sorting indices
 
354
        idxs = np.argsort(X[I])
 
355
 
 
356
        # Reorder points
 
357
        X = X[I][idxs]
 
358
 
 
359
        # Compute Vandermonde matrix
 
360
        V = P.GradVandermonde1D(X, npoint,0,norm=True)
 
361
 
 
362
        # Construct weights
 
363
        (W,res,rank,s) = npla.lstsq(V.T,b)
 
364
        # (W,res) = opt.nnls(V.T,b)
 
365
 
 
366
        # Test orthogonality of new quadrature rule
 
367
        I_tmp = np.dot(V.T,np.dot(np.diag(W),V))
 
368
        plt.figure()
 
369
        plt.imshow( np.log10(np.abs(I_tmp)), interpolation='none')
 
370
        plt.colorbar()
 
371
        plt.title("Iter %d" % it)
 
372
        plt.show(block=False)
 
373
 
 
374
        x_list.append( X[:] )
 
375
 
 
376
    # Plot Points
 
377
    plt.figure()
 
378
    for i in range(tot_it-1):
 
379
        plt.plot(x_list[i+1], i * np.ones(x_list[i+1].shape), 'go')
 
380
        plt.plot(x_list[i], i * np.ones(x_list[i].shape), 'ko')
 
381
    plt.xlim([-1.1,1.1])
 
382
    plt.ylim([-1,tot_it-1])
 
383
    plt.show(False)