Skip to content

Commit 569b9f5

Browse files
carlosmunozmoncayocarlosmunozmoncayo
committed
Working for very short times
I changed the sign of the surface term. With negative sign (or 0) in front of that term the solution blows up immediately.
1 parent c456b0e commit 569b9f5

File tree

3 files changed

+90
-48
lines changed

3 files changed

+90
-48
lines changed

‎SBP_matrices.py

Lines changed: 6 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,9 @@
11
import numpy as np
22
from basis_nodes import generate_lagrange_poly, generate_LGL_points
33

4-
def first_order_D( x_nodes='LGL' , n=10 ):
5-
if x_nodes == 'LGL':
4+
def first_order_D( x_nodes , n=10, use_LGL=False ):
5+
if use_LGL:
66
_,_,_,_,_,_,x_nodes,_ = generate_LGL_points(n-1)
7-
elif x_nodes == 'LG':
8-
_,_,_,_,x_nodes,_,_,_ = generate_LGL_points(n-1)
97

108
n_nodes = len(x_nodes)
119
D = np.zeros((n_nodes,n_nodes))
@@ -15,18 +13,14 @@ def first_order_D( x_nodes='LGL' , n=10 ):
1513
D[i,j] = Ljp(x_nodes[i])
1614
return D
1715

18-
def first_order_P_Q( x_Lagrange_nodes='LGL', x_abcissae='LGL' , n=10 ):
19-
if x_Lagrange_nodes == 'LGL':
16+
def first_order_P_Q( x_Lagrange_nodes, x_abcissae , w_abcissae, n=10, use_LGL=False ):
17+
if use_LGL:
2018
_,_,_,_,_,_,x_Lagrange_nodes,_ = generate_LGL_points(n-1)
2119

2220
n_nodes = len(x_Lagrange_nodes)
23-
24-
if x_abcissae == 'LGL':
21+
if use_LGL:
2522
_,_,_,_,_,_,x_abcissae, w_abcissae = generate_LGL_points(n_nodes-1)
26-
27-
_,_,_,_,_,_,_, w_abcissae = generate_LGL_points(n_nodes-1)
28-
w_abcissae = 0.5*(x_abcissae[-1]-x_abcissae[0])*w_abcissae
29-
23+
3024
#Generating a list of Lagrange basis
3125
list_Lagrange_poly = []
3226
list_Lagrange_poly_prime = []

‎SBP_nodes_matrices.ipynb

Lines changed: 82 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -129,8 +129,25 @@
129129
"degree_basis = 10 #degree_basis = n_points-1\n",
130130
"(p_Legn, p_Legn_prime,\n",
131131
"p_Legn1, p_Legn1_prime,\n",
132-
"x_LG, w_LG,\n",
133-
"x_LGL, w_LGL) = generate_LGL_points(degree_basis)\n",
132+
"x_LG_unsorted, w_LG_unsorted,\n",
133+
"x_LGL_unsorted, w_LGL_unsorted) = generate_LGL_points(degree_basis)\n",
134+
"\n",
135+
"#####\n",
136+
"#The root finder returns unsorted roots,\n",
137+
"#We fix that in these lines\n",
138+
"#####\n",
139+
"#Indexes that would sort x_LG using value\n",
140+
"sort_idxs_LG = np.argsort(x_LG_unsorted)\n",
141+
"#Sorting x_LG and w_LG using those indexes\n",
142+
"x_LG = x_LG_unsorted[sort_idxs_LG]\n",
143+
"w_LG = w_LG_unsorted[sort_idxs_LG]\n",
144+
"\n",
145+
"#Indexes that would sort x_LG using value\n",
146+
"sort_idxs_LGL = np.argsort(x_LGL_unsorted)\n",
147+
"#Sorting x_LG and w_LG using those indexes\n",
148+
"x_LGL = x_LGL_unsorted[sort_idxs_LGL]\n",
149+
"w_LGL = w_LGL_unsorted[sort_idxs_LGL]\n",
150+
"\n",
134151
"x = np.linspace(-1,1,1000)"
135152
]
136153
},
@@ -243,7 +260,7 @@
243260
},
244261
{
245262
"cell_type": "markdown",
246-
"id": "29b4e4b2",
263+
"id": "fd474626",
247264
"metadata": {},
248265
"source": [
249266
"### First Order Differential Operator\n",
@@ -264,8 +281,24 @@
264281
"degree_basis = 2 #degree_basis = n_nodes-1\n",
265282
"(p_Legn, p_Legn_prime,\n",
266283
"p_Legn1, p_Legn1_prime,\n",
267-
"x_LG, w_LG,\n",
268-
"x_LGL, w_LGL) = generate_LGL_points(degree_basis)\n",
284+
"x_LG_unsorted, w_LG_unsorted,\n",
285+
"x_LGL_unsorted, w_LGL_unsorted) = generate_LGL_points(degree_basis)\n",
286+
"\n",
287+
"#####\n",
288+
"#The root finder returns unsorted roots,\n",
289+
"#We fix that in these lines\n",
290+
"#####\n",
291+
"#Indexes that would sort x_LG using value\n",
292+
"sort_idxs_LG = np.argsort(x_LG_unsorted)\n",
293+
"#Sorting x_LG and w_LG using those indexes\n",
294+
"x_LG = x_LG_unsorted[sort_idxs_LG]\n",
295+
"w_LG = w_LG_unsorted[sort_idxs_LG]\n",
296+
"\n",
297+
"#Indexes that would sort x_LG using value\n",
298+
"sort_idxs_LGL = np.argsort(x_LGL_unsorted)\n",
299+
"#Sorting x_LG and w_LG using those indexes\n",
300+
"x_LGL = x_LGL_unsorted[sort_idxs_LGL]\n",
301+
"w_LGL = w_LGL_unsorted[sort_idxs_LGL]\n",
269302
"\n",
270303
"n_nodes = len(x_LGL)\n",
271304
"x = np.linspace(-0.5,0.5,n_nodes)\n",
@@ -282,7 +315,7 @@
282315
},
283316
{
284317
"cell_type": "markdown",
285-
"id": "3b0cdff7",
318+
"id": "c42c58a3",
286319
"metadata": {},
287320
"source": [
288321
"### Computing P and Q\n",
@@ -299,9 +332,9 @@
299332
"metadata": {},
300333
"outputs": [],
301334
"source": [
302-
"P_LG, Q_LG = first_order_P_Q(x_Lagrange_nodes=x_LG)\n",
303-
"P_LGL, Q_LGL = first_order_P_Q(x_Lagrange_nodes=x_LGL)\n",
304-
"P, Q = first_order_P_Q(x_Lagrange_nodes=x)\n",
335+
"P_LG, Q_LG = first_order_P_Q(x_Lagrange_nodes=x_LG, x_abcissae=x_LG, w_abcissae=w_LG)\n",
336+
"P_LGL, Q_LGL = first_order_P_Q(x_Lagrange_nodes=x_LGL, x_abcissae=x_LGL, w_abcissae=w_LGL)\n",
337+
"P, Q = first_order_P_Q(x_Lagrange_nodes=x, x_abcissae=x_LGL, w_abcissae=w_LGL)\n",
305338
"print(\"P_LG:\")\n",
306339
"print(np.round(P_LG,2))\n",
307340
"print(\"P_LGL:\")\n",
@@ -415,15 +448,15 @@
415448
},
416449
{
417450
"cell_type": "markdown",
418-
"id": "6b9af948",
451+
"id": "6b96145a",
419452
"metadata": {},
420453
"source": [
421454
"## LGL nodes and weights, and SBP operators mapped to a new interval"
422455
]
423456
},
424457
{
425458
"cell_type": "markdown",
426-
"id": "21e75e21",
459+
"id": "ba077555",
427460
"metadata": {},
428461
"source": [
429462
"It is often desired to work in a different interval than $[-1,1]$, say $[x_\\min, x_\\max]$. \n",
@@ -450,36 +483,52 @@
450483
{
451484
"cell_type": "code",
452485
"execution_count": null,
453-
"id": "5d3c6e4d",
486+
"id": "0b1296b3",
454487
"metadata": {},
455488
"outputs": [],
456489
"source": [
457490
"degree_basis = 4 #degree_basis = n_nodes-1\n",
458-
"(p_Legn, p_Legn_prime,\n",
459-
"p_Legn1, p_Legn1_prime,\n",
460-
"xi_LG, w_LG,\n",
461-
"xi_LGL, w_LGL) = generate_LGL_points(degree_basis)\n",
491+
"(_, _,\n",
492+
"_, _,\n",
493+
"_, _,\n",
494+
"xi_LGL_unsorted, w_LGL_unsorted) = generate_LGL_points(degree_basis)\n",
495+
"\n",
496+
"#######################\n",
497+
"#######################\n",
498+
"#The root finder returns unsorted roots,\n",
499+
"#We fix that in these lines\n",
500+
"\n",
501+
"#Indexes that would sort x_LG using value\n",
502+
"sort_idxs_LGL = np.argsort(xi_LGL_unsorted)\n",
503+
"#Sorting x_LG and w_LG using those indexes\n",
504+
"xi_LGL = xi_LGL_unsorted[sort_idxs_LGL]\n",
505+
"w_LGL = w_LGL_unsorted[sort_idxs_LGL]\n",
462506
"\n",
463507
"xmin = -3; xmax=3\n",
464508
"a = (xmax-xmin)/2.; b = (xmax+xmin)/2.\n",
465509
"x_mapped = a*xi_LGL+b\n",
466-
"w_mapped = a*w_LGL"
510+
"w_mapped = a*w_LGL\n",
511+
"#######################\n",
512+
"#######################\n",
513+
"\n"
467514
]
468515
},
469516
{
470517
"cell_type": "code",
471518
"execution_count": null,
472-
"id": "eab7e897",
519+
"id": "89b8b473",
473520
"metadata": {},
474521
"outputs": [],
475522
"source": [
523+
"print(f\"x={x_mapped}\")\n",
524+
"print(f\"Sum of weights\")\n",
476525
"print(np.sum(w_LGL))\n",
477526
"print(np.sum(w_mapped))"
478527
]
479528
},
480529
{
481530
"cell_type": "markdown",
482-
"id": "d5880876",
531+
"id": "51b660f3",
483532
"metadata": {},
484533
"source": [
485534
"Since the Lagrange basis and, therefore, the first order differential operator just require a set of nodes for computation, they are computed exactly as before.\n",
@@ -490,28 +539,29 @@
490539
{
491540
"cell_type": "code",
492541
"execution_count": null,
493-
"id": "6b5647d5",
542+
"id": "f9eb7843",
494543
"metadata": {},
495544
"outputs": [],
496545
"source": [
497-
"poly, polyp = generate_lagrange_poly(1,x_mapped)\n",
546+
"idx_support = 3 #Which node is going to be 1\n",
547+
"poly, polyp = generate_lagrange_poly(idx_support,x_mapped)\n",
498548
"x = np.linspace(-3.5,3.5,50)\n",
499549
"vpoly = np.vectorize(poly)\n",
500550
"vpolyp = np.vectorize(polyp)\n",
501551
"fig = plt.figure(dpi=100, figsize=(8,4))\n",
502552
"ax = fig.add_subplot(1, 1, 1)\n",
503553
"ax.set_ylim([-5,5])\n",
504554
"ax.plot(x,np.zeros(len(x)),color='k', lw=1)\n",
505-
"ax.plot(x, vpoly(x), label=f'$L_{1}(x)$')\n",
506-
"ax.plot(x, vpolyp(x),color='g', label=f'$L_{1}\\'(x)$')\n",
555+
"ax.plot(x, vpoly(x), label=f'$L_{idx_support}(x)$')\n",
556+
"ax.plot(x, vpolyp(x),color='g', label=f'$L_{idx_support}\\'(x)$')\n",
507557
"ax.scatter(x_mapped,vpoly(x_mapped),marker='x', color='r')\n",
508558
"plt.legend()\n",
509559
"plt.show()"
510560
]
511561
},
512562
{
513563
"cell_type": "markdown",
514-
"id": "8f096353",
564+
"id": "b8f02bdd",
515565
"metadata": {},
516566
"source": [
517567
"### First Order Differential Operator D"
@@ -520,7 +570,7 @@
520570
{
521571
"cell_type": "code",
522572
"execution_count": null,
523-
"id": "8a818d8b",
573+
"id": "4507dee6",
524574
"metadata": {},
525575
"outputs": [],
526576
"source": [
@@ -530,7 +580,7 @@
530580
},
531581
{
532582
"cell_type": "markdown",
533-
"id": "b7d79ad0",
583+
"id": "a8b92631",
534584
"metadata": {},
535585
"source": [
536586
"### Computing P and Q\n",
@@ -545,13 +595,13 @@
545595
{
546596
"cell_type": "code",
547597
"execution_count": null,
548-
"id": "4e4dde46",
598+
"id": "9f690af2",
549599
"metadata": {
550600
"scrolled": true
551601
},
552602
"outputs": [],
553603
"source": [
554-
"P, Q = first_order_P_Q(x_Lagrange_nodes=x_mapped, x_abcissae=x_mapped)\n",
604+
"P, Q = first_order_P_Q(x_Lagrange_nodes=x_mapped, x_abcissae=x_mapped, w_abcissae=w_mapped)\n",
555605
"print(\"P:\")\n",
556606
"print(np.round(P,2))\n",
557607
"print(\"Q:\")\n",
@@ -560,7 +610,7 @@
560610
},
561611
{
562612
"cell_type": "markdown",
563-
"id": "51c85785",
613+
"id": "585fcfa2",
564614
"metadata": {},
565615
"source": [
566616
"### Verifying SBP properties\n",
@@ -570,7 +620,7 @@
570620
{
571621
"cell_type": "code",
572622
"execution_count": null,
573-
"id": "ce1d4970",
623+
"id": "c3b4dc89",
574624
"metadata": {},
575625
"outputs": [],
576626
"source": [
@@ -579,7 +629,7 @@
579629
},
580630
{
581631
"cell_type": "markdown",
582-
"id": "33b9a684",
632+
"id": "1ccd0d8b",
583633
"metadata": {},
584634
"source": [
585635
"Now we check if\n",
@@ -591,7 +641,7 @@
591641
{
592642
"cell_type": "code",
593643
"execution_count": null,
594-
"id": "f6ad6b75",
644+
"id": "46f1b6a5",
595645
"metadata": {},
596646
"outputs": [],
597647
"source": [

‎basis_nodes.py

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,9 @@
99
#print(cubic(x))
1010
#print(2.5*x**3-1.5*x)
1111

12-
def generate_lagrange_poly(j, x_nodes='LGL',n=10):
13-
if x_nodes == 'LGL':
12+
def generate_lagrange_poly(j, x_nodes,n=10, use_LGL=False):
13+
if use_LGL:
1414
_,_,_,_,_,_,x_nodes,_= generate_LGL_points(n-1)
15-
elif x_nodes == 'LG':
16-
_,_,_,_,x_nodes,_,_,_= generate_LGL_points(n-1)
1715
xj = x_nodes[j]
1816
n_nodes = len(x_nodes)
1917
def Lj(x):

0 commit comments

Comments
 (0)