hedge-notes.tm 82.83 KiB
<TeXmacs|1.0.7>
<style|<tuple|article|maxima|axiom|mystyle>>
<\body>
<section|Mapping from 2D equilateral to Unit Coordinates>
<with|prog-language|axiom|prog-session|default|<\session>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
)clear all
</input>
<\output>
\ \ \ All user variables and function definitions have been cleared.
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
verteq:=[2/sqrt(3)*vector [sin(i*2*%pi/3),cos(i*2*%pi/3)] for i in
0..2]
</input>
<\output>
<with|mode|math|math-display|true|<left|[><left|[>0,<space|0.5spc><frac|2<sqrt|3>|3><right|]>,<space|0.5spc><left|[>1,<space|0.5spc>-<frac|<sqrt|3>|3><right|]>,<space|0.5spc><left|[>-1,<space|0.5spc>-<frac|<sqrt|3>|3><right|]><right|]><leqno>(1)>
<axiomtype|List Vector Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
vertun:= [vector[-1,1], vector[1,-1], vector[-1,-1]]
</input>
<\output>
<with|mode|math|math-display|true|<left|[><left|[>-1,<space|0.5spc>1<right|]>,<space|0.5spc><left|[>1,<space|0.5spc>-1<right|]>,<space|0.5spc><left|[>-1,<space|0.5spc>-1<right|]><right|]><leqno>(2)>
<axiomtype|List Vector Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
A:=matrix [[a,b],[c,d]]
</input>
<\output>
<with|mode|math|math-display|true|<left|[><tabular*|<tformat|<cwith|1|-1|1|1|cell-halign|c>|<cwith|1|-1|2|2|cell-halign|c>|<table|<row|<cell|a>|<cell|b>>|<row|<cell|c>|<cell|d>>>>><right|]><leqno>(3)>
<axiomtype|Matrix Polynomial Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
v:=vector [e,f]
</input>
<\output>
<with|mode|math|math-display|true|<left|[>e,<space|0.5spc>f<right|]><leqno>(4)>
<axiomtype|Vector OrderedVariableList [e,f] >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
eqns:=concat(
[(A*verteq.i+v).1=vertun.i.1 for i in 1..3],
[(A*verteq.i+v).2=vertun.i.2 for i in 1..3])
</input>
<\output>
<with|mode|math|math-display|true|<left|[><frac|2b<sqrt|3>+3e|3>=-1,<space|0.5spc><frac|-b<sqrt|3>+3e+3a|3>=1,<space|0.5spc><frac|-b<sqrt|3>+3e-3a|3>=-1,<space|0.5spc><frac|2d<sqrt|3>+3f|3>=1,<space|0.5spc><frac|-d<sqrt|3>+3f+3c|3>=-1,<space|0.5spc><frac|-d<sqrt|3>+3f-3c|3>=-1<right|]><leqno>(5)>
<axiomtype|List Equation Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
solve(eqns,[a,b,c,d,e,f])
</input>
<\output>
<with|mode|math|math-display|true|<left|[><left|[>a=1,<space|0.5spc>b=-<frac|1|<sqrt|3>>,<space|0.5spc>c=0,<space|0.5spc>d=<frac|2|<sqrt|3>>,<space|0.5spc>e=-<frac|1|3>,<space|0.5spc>f=-<frac|1|3><right|]><right|]><leqno>(6)>
<axiomtype|List List Equation Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
\;
</input>
</session>>
<section|Mapping from 3D equilateral to Unit Coordinates>
<with|prog-language|axiom|prog-session|default|<\session>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
)clear all
</input>
<\output>
\ \ \ All user variables and function definitions have been cleared.
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
verteq:=[
vector [-1,-1/sqrt(3),-1/sqrt(6)],
vector [ 1,-1/sqrt(3),-1/sqrt(6)],
vector [ 0, 2/sqrt(3),-1/sqrt(6)],
vector [ 0, \ \ \ \ \ \ \ \ 0, 3/sqrt(6)]
]
</input>
<\output>
<with|mode|math|math-display|true|<left|[><left|[>-1,<space|0.5spc>-<frac|<sqrt|3>|3>,<space|0.5spc>-<frac|<sqrt|6>|6><right|]>,<space|0.5spc><left|[>1,<space|0.5spc>-<frac|<sqrt|3>|3>,<space|0.5spc>-<frac|<sqrt|6>|6><right|]>,<space|0.5spc><left|[>0,<space|0.5spc><frac|2<sqrt|3>|3>,<space|0.5spc>-<frac|<sqrt|6>|6><right|]>,<space|0.5spc><left|[>0,<space|0.5spc>0,<space|0.5spc><frac|<sqrt|6>|2><right|]><right|]><leqno>(1)>
<axiomtype|List Vector AlgebraicNumber >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
vertun:= [
[-1,-1,-1],
[ 1,-1,-1],
[-1, 1,-1],
[-1,-1, 1]
]
</input>
<\output>
<with|mode|math|math-display|true|<left|[><left|[>-1,<space|0.5spc>-1,<space|0.5spc>-1<right|]>,<space|0.5spc><left|[>1,<space|0.5spc>-1,<space|0.5spc>-1<right|]>,<space|0.5spc><left|[>-1,<space|0.5spc>1,<space|0.5spc>-1<right|]>,<space|0.5spc><left|[>-1,<space|0.5spc>-1,<space|0.5spc>1<right|]><right|]><leqno>(2)>
<axiomtype|List List Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
A:=matrix [[a,b,c],[d,e,f],[g,h,i]]
</input>
<\output>
<with|mode|math|math-display|true|<left|[><tabular*|<tformat|<cwith|1|-1|1|1|cell-halign|c>|<cwith|1|-1|2|2|cell-halign|c>|<cwith|1|-1|3|3|cell-halign|c>|<table|<row|<cell|a>|<cell|b>|<cell|c>>|<row|<cell|d>|<cell|e>|<cell|f>>|<row|<cell|g>|<cell|h>|<cell|i>>>>><right|]><leqno>(3)>
<axiomtype|Matrix Polynomial Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
v:=vector [j,k,l]
</input>
<\output>
<with|mode|math|math-display|true|<left|[>j,<space|0.5spc>k,<space|0.5spc>l<right|]><leqno>(4)>
<axiomtype|Vector OrderedVariableList [j,k,l] >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
eqns:=concat(
[(A*verteq.i+v).1=vertun.i.1 for i in 1..4],
concat(
[(A*verteq.i+v).2=vertun.i.2 for i in 1..4],
[(A*verteq.i+v).3=vertun.i.3 for i in 1..4]
))
</input>
<\output>
<with|mode|math|math-display|true|<left|[>j-<frac|<sqrt|6>|6>c-<frac|<sqrt|3>|3>b-a=-1,<space|0.5spc>j-<frac|<sqrt|6>|6>c-<frac|<sqrt|3>|3>b+a=1,<space|0.5spc>j-<frac|<sqrt|6>|6>c+<frac|2<sqrt|3>|3>b=-1,<space|0.5spc>j+<frac|<sqrt|6>|2>c=-1,<space|0.5spc>k-<frac|<sqrt|6>|6>f-<frac|<sqrt|3>|3>e-d=-1,<space|0.5spc>k-<frac|<sqrt|6>|6>f-<frac|<sqrt|3>|3>e+d=-1,<space|0.5spc>k-<frac|<sqrt|6>|6>f+<frac|2<sqrt|3>|3>e=1,<space|0.5spc>k+<frac|<sqrt|6>|2>f=-1,<space|0.5spc>l-<frac|<sqrt|6>|6>i-<frac|<sqrt|3>|3>h-g=-1,<space|0.5spc>l-<frac|<sqrt|6>|6>i-<frac|<sqrt|3>|3>h+g=-1,<space|0.5spc>l-<frac|<sqrt|6>|6>i+<frac|2<sqrt|3>|3>h=-1,<space|0.5spc>l+<frac|<sqrt|6>|2>i=1<right|]><leqno>(5)>
<axiomtype|List Equation Polynomial AlgebraicNumber >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
solve(eqns,[a,b,c,d,e,f,g,h,i,j,k,l])
</input>
<\output>
<with|mode|math|math-display|true|<left|[><left|[>a=1,<space|0.5spc>b=-<frac|<sqrt|3>|3>,<space|0.5spc>c=-<frac|<sqrt|6>|6>,<space|0.5spc>d=0,<space|0.5spc>e=<frac|2<sqrt|3>|3>,<space|0.5spc>f=-<frac|<sqrt|6>|6>,<space|0.5spc>g=0,<space|0.5spc>h=0,<space|0.5spc>i=<frac|<sqrt|6>|2>,<space|0.5spc>j=-<frac|1|2>,<space|0.5spc>k=-<frac|1|2>,<space|0.5spc>l=-<frac|1|2><right|]><right|]><leqno>(6)>
<axiomtype|List List Equation Fraction Polynomial AlgebraicNumber >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
\;
</input>
</session>>
<section|Derivatives of the 2D Basis functions>
<with|prog-language|axiom|prog-session|default|<\session>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
)clear all
</input>
<\output>
\ \ \ All user variables and function definitions have been cleared.
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
f:=operator 'f
</input>
<\output>
<with|mode|math|math-display|true|f<leqno>(1)>
<axiomtype|BasicOperator >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
g:=operator 'g
</input>
<\output>
<with|mode|math|math-display|true|g<leqno>(2)>
<axiomtype|BasicOperator >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
a:=2*(1+r)/(1-s)-1
</input>
<\output>
<with|mode|math|math-display|true|<frac|-s-2r-1|s-1><leqno>(3)>
<axiomtype|Fraction Polynomial Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
D(a,r)
</input>
<\output>
<with|mode|math|math-display|true|-<frac|2|s-1><leqno>(4)>
<axiomtype|Fraction Polynomial Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
D(a,s)
</input>
<\output>
<with|mode|math|math-display|true|<frac|2r+2|s<rsup|2>-2s+1><leqno>(5)>
<axiomtype|Fraction Polynomial Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
expr:=sqrt(2)*f(a)*g(s)*(1-s)^i
</input>
<\output>
<with|mode|math|math-display|true|<sqrt|2>f<left|(><frac|-s-2r-1|s-1><right|)>g<left|(>s<right|)><left|(>-s+1<right|)><rsup|i><leqno>(6)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
D(expr,r)
</input>
<\output>
<with|mode|math|math-display|true|-<frac|2<sqrt|2>g<left|(>s<right|)><left|(>-s+1<right|)><rsup|i>f<rsub|
><rsup|,><left|(><frac|-s-2r-1|s-1><right|)>|s-1><leqno>(7)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
D(expr,s)
*
(s-1)^2 </input>
<\output>
<with|mode|math|math-display|true|<left|(>s<rsup|2>-2s+1<right|)><sqrt|2>f<left|(><frac|-s-2r-1|s-1><right|)><left|(>-s+1<right|)><rsup|i>g<rsub|
><rsup|,><left|(>s<right|)>+<left|(>2r+2<right|)><sqrt|2>g<left|(>s<right|)><left|(>-s+1<right|)><rsup|i>f<rsub|
><rsup|,><left|(><frac|-s-2r-1|s-1><right|)>+<left|(>-i*s<rsup|2>+2i*s-i<right|)><sqrt|2>f<left|(><frac|-s-2r-1|s-1><right|)>g<left|(>s<right|)><left|(>-s+1<right|)><rsup|<left|(>i-1<right|)>><leqno>(9)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
\;
</input>
</session>>
<\eqnarray*>
<tformat|<table|<row|<cell|>|<cell|>|<cell|<frac|1|(s-1)<rsup|2>><left|[><left|(>s<rsup|2>-2s+1<right|)><sqrt|2>f<left|(><frac|-s-2r-1|s-1><right|)><left|(>-s+1<right|)><rsup|i>g<rsub|
><rsup|,><left|(>s<right|)><right|]>>>|<row|<cell|>|<cell|>|<cell|+<frac|1|(s-1)<rsup|2>><left|[><left|(>2r+2<right|)><sqrt|2>g<left|(>s<right|)><left|(>-s+1<right|)><rsup|i>f<rsub|
><rsup|,><left|(><frac|-s-2r-1|s-1><right|)><right|]>>>|<row|<cell|>|<cell|>|<cell|+<frac|1|(s-1)<rsup|2>><left|[><left|(>-i*s<rsup|2>+2i*s-i<right|)><sqrt|2>f<left|(><frac|-s-2r-1|s-1><right|)>g<left|(>s<right|)><left|(>-s+1<right|)><rsup|<left|(>i-1<right|)>><right|]>>>|<row|<cell|>|<cell|=>|<cell|<sqrt|2><left|[>f<left|(>a<right|)><left|(>1-s<right|)><rsup|i>g<rsub|
><rsup|,><left|(>s<right|)>>>|<row|<cell|>|<cell|>|<cell|+<left|(>2r+2<right|)>g<left|(>s<right|)><left|(>1-s<right|)><rsup|i-2>f<rsub|
><rsup|,><left|(>a<right|)>>>|<row|<cell|>|<cell|>|<cell|-i*f<left|(>a<right|)>g<left|(>s<right|)><left|(>1-s<right|)><rsup|<left|(>i-1<right|)>><right|]>>>>>
</eqnarray*>
<section|Derivatives of the 3D Basis functions>
<with|prog-language|axiom|prog-session|default|<\session>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
)clear all
</input>
<\output>
\ \ \ All user variables and function definitions have been cleared.
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
f:=operator 'f
</input>
<\output>
<with|mode|math|math-display|true|f<leqno>(1)>
<axiomtype|BasicOperator >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
g:=operator 'g
</input>
<\output>
<with|mode|math|math-display|true|g<leqno>(2)>
<axiomtype|BasicOperator >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
h:=operator 'h
</input>
<\output>
<with|mode|math|math-display|true|h<leqno>(3)>
<axiomtype|BasicOperator >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
a:= -2*(1+r)/(s+t) - 1
</input>
<\output>
<with|mode|math|math-display|true|<frac|-t-s-2r-2|t+s><leqno>(4)>
<axiomtype|Fraction Polynomial Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
b:=2*(1+s)/(1-t) - 1
</input>
<\output>
<with|mode|math|math-display|true|<frac|-t-2s-1|t-1><leqno>(5)>
<axiomtype|Fraction Polynomial Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
expr:=sqrt(8)*f(a)*g(b)*(1-b)^i*h(c)*(1-c)**(i+j)
</input>
<\output>
<with|mode|math|math-display|true|2<sqrt|2>f<left|(><frac|-t-s-2r-2|t+s><right|)>g<left|(><frac|-t-2s-1|t-1><right|)>h<left|(>c<right|)><left|(>-c+1<right|)><rsup|<left|(>j+i<right|)>><frac|2t+2s|t-1><rsup|i><leqno>(6)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
D(expr,r)
</input>
<\output>
<with|mode|math|math-display|true|-<frac|4<sqrt|2>g<left|(><frac|-t-2s-1|t-1><right|)>h<left|(>c<right|)><left|(>-c+1<right|)><rsup|<left|(>j+i<right|)>><frac|2t+2s|t-1><rsup|i>f<rsub|
><rsup|,><left|(><frac|-t-s-2r-2|t+s><right|)>|t+s><leqno>(7)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
dfds:=D(expr,s);
</input>
<\output>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
dfdsden:=denominator dfds
</input>
<\output>
<with|mode|math|math-display|true|t<rsup|3>+<left|(>2s-1<right|)>t<rsup|2>+<left|(>s<rsup|2>-2s<right|)>t-s<rsup|2><leqno>(20)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
dfds*dfdsden
</input>
<\output>
<with|mode|math|math-display|true|<left|(>-4t<rsup|2>-8s*t-4s<rsup|2><right|)><sqrt|2>f<left|(><frac|-t-s-2r-2|t+s><right|)>h<left|(>c<right|)><left|(>-c+1<right|)><rsup|<left|(>j+i<right|)>><frac|2t+2s|t-1><rsup|i>g<rsub|
><rsup|,><left|(><frac|-t-2s-1|t-1><right|)>+<left|(><left|(>4r+4<right|)>t-4r-4<right|)><sqrt|2>g<left|(><frac|-t-2s-1|t-1><right|)>h<left|(>c<right|)><left|(>-c+1<right|)><rsup|<left|(>j+i<right|)>><frac|2t+2s|t-1><rsup|i>f<rsub|
><rsup|,><left|(><frac|-t-s-2r-2|t+s><right|)>+<left|(>4i*t<rsup|2>+8i*s*t+4i*s<rsup|2><right|)><sqrt|2>f<left|(><frac|-t-s-2r-2|t+s><right|)>g<left|(><frac|-t-2s-1|t-1><right|)>h<left|(>c<right|)><left|(>-c+1<right|)><rsup|<left|(>j+i<right|)>><frac|2t+2s|t-1><rsup|<left|(>i-1<right|)>><leqno>(18)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
dfdt:=D(expr,t);
</input>
<\output> <axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
dfdtden:=denominator dfdt
</input>
<\output>
<with|mode|math|math-display|true|t<rsup|4>+<left|(>2s-2<right|)>t<rsup|3>+<left|(>s<rsup|2>-4s+1<right|)>t<rsup|2>+<left|(>-2s<rsup|2>+2s<right|)>t+s<rsup|2><leqno>(23)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
dfdt*dfdtden
</input>
<\output>
<with|mode|math|math-display|true|<left|(><left|(>4s+4<right|)>t<rsup|2>+<left|(>8s<rsup|2>+8s<right|)>t+4s<rsup|3>+4s<rsup|2><right|)><sqrt|2>f<left|(><frac|-t-s-2r-2|t+s><right|)>h<left|(>c<right|)><left|(>-c+1<right|)><rsup|<left|(>j+i<right|)>><frac|2t+2s|t-1><rsup|i>g<rsub|
><rsup|,><left|(><frac|-t-2s-1|t-1><right|)>+<left|(><left|(>4r+4<right|)>t<rsup|2>+<left|(>-8r-8<right|)>t+4r+4<right|)><sqrt|2>g<left|(><frac|-t-2s-1|t-1><right|)>h<left|(>c<right|)><left|(>-c+1<right|)><rsup|<left|(>j+i<right|)>><frac|2t+2s|t-1><rsup|i>f<rsub|
><rsup|,><left|(><frac|-t-s-2r-2|t+s><right|)>+<left|(><left|(>-4i*s-4i<right|)>t<rsup|2>+<left|(>-8i*s<rsup|2>-8i*s<right|)>t-4i*s<rsup|3>-4i*s<rsup|2><right|)><sqrt|2>f<left|(><frac|-t-s-2r-2|t+s><right|)>g<left|(><frac|-t-2s-1|t-1><right|)>h<left|(>c<right|)><left|(>-c+1<right|)><rsup|<left|(>j+i<right|)>><frac|2t+2s|t-1><rsup|<left|(>i-1<right|)>><leqno>(24)>
<axiomtype|Expression Integer >
</output>
<\input|<with|color|red|<with|mode|math|\<rightarrow\>> >>
\;
</input>
</session>>
<section|DG Schemes>
<subsection|Notation>
Let <with|mode|math|l<rsub|i>> be the <with|mode|math|i>th Lagrange
interpolation polynomial and <with|mode|math|\<b-u\>=(u<rsub|i>)> be the
vector of nodal coefficients of our approximated function
<with|mode|math|<wide|u|~>(x)=u<rsub|i>l<rsub|i>(x)>. (<with|mode|math|u>
without subscript denotes a continuous function <with|mode|math|u>.)
<with|mode|math|F\<subset\>\<partial\>T<rsub|k>> represents the discrete
faces of the element <with|mode|math|T<rsub|k>>, and
<with|mode|math|l<rsub|i><rsup|F>> is the Lagrange interpolation polynomial
on the face <with|mode|math|f> corresponding to the node with global number
<with|mode|math|i>, or constant zero if there is no such function (such as
when <with|mode|math|i> is not the number of a node in <with|mode|math|F>).
Further,
<\eqnarray*>
<tformat|<table|<row|<cell|M<rsub|i j><rsup|k>>|<cell|\<assign\>>|<cell|<big|int><rsub|T<rsub|k>>l<rsub|i>l<rsub|j>,>>|<row|<cell|S<rsub|i
j,\<partial\>x<rsub|\<nu\>>><rsup|k>>|<cell|\<assign\>>|<cell|<big|int><rsub|T<rsub|k>>l<rsub|i>\<partial\><rsub|x<rsub|\<nu\>>>l<rsub|j>,>>|<row|<cell|M<rsub|i
j><rsup|F>>|<cell|\<assign\>>|<cell|<big|int><rsub|F>l<rsub|i><rsup|F>l<rsub|j><rsup|F>.>>>>
</eqnarray*>
Note that <with|mode|math|M<rsup|k>=(M<rsup|k>)<rsup|T>> and
<with|mode|math|M<rsup|F>=(M<rsup|F>)<rsup|T>>.
Recall <with|mode|math|V<wide|\<b-u\>|^>=\<b-u\>> where
<with|mode|math|<wide|u|^><rsub|i>p<rsub|i>=u<rsub|i>l<rsub|i>>, and
<\equation*>
V<rsup|T>\<b-l\>(x)=<matrix|<tformat|<table|<row|<cell|p<rsub|1>(x<rsub|1>)>|<cell|>|<cell|p(x<rsub|n>)>>|<row|<cell|\<vdots\>>|<cell|>|<cell|\<vdots\>>>|<row|<cell|p<rsub|n>(x<rsub|1>)>|<cell|\<cdots\>>|<cell|p<rsub|n>(x<rsub|n>)>>>>>\<b-l\>(x)=<matrix|<tformat|<table|<row|<cell|l<rsub|1>(x)p<rsub|1>(x<rsub|1>)+\<cdots\>+l<rsub|n>(x)p<rsub|1>(x<rsub|n>)>>|<row|<cell|\<vdots\>>>|<row|<cell|l<rsub|1>(x)p<rsub|n>(x<rsub|1>)+\<cdots\>+l<rsub|n>(x)p<rsub|n>(x<rsub|n>)>>>>>=<matrix|<tformat|<table|<row|<cell|p<rsub|1>(x)>>|<row|<cell|\<vdots\>>>|<row|<cell|p<rsub|n>(x)>>>>>,
</equation*>
simply because the sum of Lagrange polynomials uniquely interpolates
<with|mode|math|p<rsub|i>(x)>. We thus find
<\eqnarray*>
<tformat|<table|<row|<cell|M<rsup|k><rsub|i
j>>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>l<rsub|i>l<rsub|j>=J<rsub|k><big|int><rsub|T>([V<rsup|-T>]<rsub|i m>p<rsub|m>)([V<rsup|-T>]<rsub|j n>p<rsub|n>)>>|<row|<cell|>|<cell|=>|<cell|J<rsub|k>[V<rsup|-T>]<rsub|i
m>[V<rsup|-T>]<rsub|j n><big|int><rsub|T>p<rsub|m>p<rsub|n>>>|<row|<cell|>|<cell|=>|<cell|J<rsub|k>[V<rsup|-T>]<rsub|i
m>[V<rsup|-T>]<rsub|j n>\<delta\><rsub|m,n>>>|<row|<cell|>|<cell|=>|<cell|J<rsub|k>[V<rsup|-T>]<rsub|i
n>[V<rsup|-T>]<rsub|j n>>>|<row|<cell|>|<cell|=>|<cell|J<rsub|k>[V<rsup|-T>]<rsub|i
n>[V<rsup|-1>]<rsub|n j>>>|<row|<cell|>|<cell|=>|<cell|J<rsub|k>[(V*V<rsup|T>)<rsup|-1>]<rsub|i
j>.>>>>
</eqnarray*>
To obtain the stiffness matrix <with|mode|math|S<rsub|i
j,\<partial\>x<rsub|\<nu\>>>>, realize that we really only need a
differentiation matrix <with|mode|math|D<rsub|x<rsub|\<nu\>>>> that maps
the nodal values of a function to those of its derivative. If we have
<with|mode|math|D<rsub|x<rsub|\<nu\>>>>, we can simply recycle the mass
matrix:
<\equation*>
S<rsub|i,j,\<partial\>x<rsub|\<nu\>>>=M*D<rsub|x<rsub|\<nu\>>>.
</equation*>
Before finding the global differentiation matrix
<with|mode|math|D<rsub|x<rsub|\<nu\>>>>, we first endeavor to find the
local differentiation matrix <with|mode|math|D<rsub|r<rsub|\<nu\>>>>, where
we denote local (unit) coordinates <with|mode|math|r> and global
coordinates <with|mode|math|x>. We define the derivative Vandermonde
matrices <with|mode|math|V<rsub|\<partial\>r<rsub|\<nu\>>>\<assign\>[\<partial\><rsub|\<nu\>>p<rsub|j>(r<rsub|i>)]<rsub|i
j><rsub|>>. Recalling <with|mode|math|V<wide|\<b-u\>|^>=\<b-u\>>, we obtain
<with|mode|math|V<rsub|\<partial\>r<rsub|\<nu\>>><wide|\<b-u\>|^>=V<rsub|\<partial\>r<rsub|\<nu\>>>V<rsup|-1>\<b-u\>=D<rsub|r<rsub|\<nu\>>>\<b-u\>>
and thus <with|mode|math|D<rsub|r<rsub|\<nu\>>>=V<rsub|\<partial\>r<rsub|\<nu\>>>V<rsup|-1>>.
For clarity, and again using the summation convention, consider the
identity
<\equation*>
<wide|u|~><rprime|'>(x)=u<rsub|j>l<rsub|j><rprime|'>(x)=u<rsub|j><wide*|l<rsub|j><rprime|'>(x<rsup|(i)>)|\<wide-underbrace\>><rsub|D<rsub|i,j>>l<rsub|i>(x).
</equation*>
If we define a function <with|mode|math|F<rsub|k>:r\<mapsto\>x>, then the
global differentiation matrix is then given by
<\equation*>
[D<rsub|x<rsub|\<nu\>>>]<rsub|i,j>=\<partial\><rsub|x<rsub|\<nu\>>>(l<rsub|j>(F<rsup|-1><rsub|k>(x<rsup|(i)>)))=(\<partial\><rsub|r>l<rsub|j>)(r<rsup|(i)><rsup|>)<left|[>(F<rsub|k><rsup|-1>)<rprime|'>(x<rsup|(i)>)<right|]><rsub|\<mu\>,\<nu\>>=[D<rsub|r<rsub|\<mu\>>>]<rsub|i,j><left|[>(F<rsup|-1><rsub|k>)<rprime|'>(x<rsup|(i)>)<right|]><rsub|\<mu\>,\<nu\>>.
</equation*>
<subsection|Advection Equation>
<subsubsection|Weak DG>
We're discretizing <with|mode|math|u<rsub|t>+v\<cdot\>\<nabla\>u=0>. The
analytic solution is
<\equation*>
u(x,t)=u<rsub|0>(x-v*t).
</equation*>
Recall
<\eqnarray*>
<tformat|<table|<row|<cell|\<nabla\>\<cdot\>(v*u*\<varphi\>)>|<cell|=>|<cell|\<nabla\>\<cdot\>(v*u)\<varphi\>+v*u\<cdot\>\<nabla\>\<varphi\>=(v\<cdot\>\<nabla\>*u)\<varphi\>+v*u\<cdot\>\<nabla\>\<varphi\>.>>>>
</eqnarray*>
Using the summation convention, we find
<\eqnarray*>
<tformat|<table|<row|<cell|0>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|T<rsub|k>>(v\<cdot\>\<nabla\>u)\<varphi\>>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>v*u\<cdot\>\<nabla\>\<varphi\>+<big|int><rsub|T<rsub|k>>\<nabla\>\<cdot\>(v*u*\<varphi\>)>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>v*u\<cdot\>\<nabla\>\<varphi\>+<big|int><rsub|\<partial\>T<rsub|k>>v*u\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|\<approx\>>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>v*u\<cdot\>\<nabla\>\<varphi\>+<big|int><rsub|\<partial\>T<rsub|k>>(v*u)<rsup|\<ast\>>\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|\<approx\>>|<cell|<big|int><rsub|T<rsub|k>>\<partial\><rsub|t>u<rsub|i>l<rsub|i>l<rsub|j>-<big|int><rsub|T<rsub|k>><matrix|<tformat|<table|<row|<cell|a<rsub|1>u<rsub|i>l<rsub|i>>>|<row|<cell|a<rsub|2>u<rsub|i>l<rsub|i>>>>>>\<cdot\><matrix|<tformat|<table|<row|<cell|\<partial\><rsub|x<rsub|2>>l<rsub|j>>>|<row|<cell|\<partial\><rsub|x<rsub|2>>l<rsub|j>>>>>>+<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><big|int><rsub|F><matrix|<tformat|<table|<row|<cell|(v<rsub|1>u<rsub|i>)<rsup|\<ast\>>l<rsub|i><rsup|F>l<rsub|j><rsup|F>>>|<row|<cell|(v<rsub|2>u<rsub|i>)<rsup|\<ast\>>l<rsub|i><rsup|F>l<rsub|j><rsup|F>>>>>>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|>>|<row|<cell|>|<cell|=>|<cell|>>>>
</eqnarray*>
<subsubsection|Strong DG>
<\eqnarray*>
<tformat|<table|<row|<cell|0>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>a*u\<cdot\>\<nabla\>\<varphi\>+<big|int><rsub|\<partial\>T<rsub|k>>(v*u)<rsup|\<ast\>>\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|D<rsub|k>>(v\<cdot\>\<nabla\>u)\<varphi\>-<big|int><rsub|\<partial\>T<rsub|k>>(v*u-{v*u})\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|D<rsub|k>>(v\<cdot\>\<nabla\>u)\<varphi\>-<big|int><rsub|\<partial\>T<rsub|k>>v<matrix|<tformat|<table|<row|<cell|<frac|1|2>>>|<row|<cell|-<frac|1|2>>>>>>\<cdot\><matrix|<tformat|<table|<row|<cell|u<rsup|->>>|<row|<cell|u<rsup|+>>>>>>\<varphi\>\<cdot\>n>>>>
</eqnarray*>
<subsection|Wave Equation>
<subsubsection|Weak DG>
\;
We have <with|mode|math|u<rsub|t t>-\<Delta\>u=0>, so
<\eqnarray*>
<tformat|<table|<row|<cell|u<rsub|t>-\<nabla\>\<cdot\>v>|<cell|=>|<cell|0,>>|<row|<cell|v<rsub|t>-\<nabla\>u>|<cell|=>|<cell|0.>>>>
</eqnarray*>
Note that <with|mode|math|v> is a vector. \ Then, using the summation
convention,
<\eqnarray*>
<tformat|<table|<row|<cell|0>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|T<rsub|k>>(\<nabla\>\<cdot\>v)\<varphi\>>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|D<rsub|k>>v\<cdot\>\<nabla\>\<varphi\>-<big|int><rsub|T<rsub|k>>\<nabla\>\<cdot\>(v\<varphi\>)>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|D<rsub|k>>v\<cdot\>\<nabla\>\<varphi\>-<big|int><rsub|\<partial\>T<rsub|k>>v\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|\<approx\>>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|D<rsub|k>>v\<cdot\>\<nabla\>\<varphi\>-<big|int><rsub|\<partial\>T<rsub|k>>v<rsup|\<ast\>>\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|\<approx\>>|<cell|<big|int><rsub|T<rsub|k>>\<partial\><rsub|t>u<rsub|i>l<rsub|i>l<rsub|j>+<big|int><rsub|T<rsub|k>><matrix|<tformat|<table|<row|<cell|[v<rsub|i>]<rsub|1>l<rsub|i>>>|<row|<cell|[v<rsub|i>]<rsub|2>l<rsub|i>>>>>>\<cdot\><matrix|<tformat|<table|<row|<cell|\<partial\><rsub|x<rsub|2>>l<rsub|j>>>|<row|<cell|\<partial\><rsub|x<rsub|2>>l<rsub|j>>>>>>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><big|int><rsub|F><matrix|<tformat|<table|<row|<cell|[v<rsub|i>]<rsup|\<ast\>><rsub|1>l<rsub|i><rsup|F>l<rsub|j><rsup|F>>>|<row|<cell|[v<rsub|i>]<rsup|\<ast\>><rsub|2>l<rsub|i><rsup|F>l<rsub|i><rsup|F>>>>>>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|M<rsub|i
j><rsup|k>\<partial\><rsub|t>u<rsub|i>+S<rsub|i
j,\<partial\>x<rsub|1>><rsup|k>[v<rsub|i>]<rsub|1>+S<rsub|i
j,\<partial\>x<rsub|2>><rsup|k>[v<rsub|i>]<rsub|2>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><left|[>[v<rsub|i>]<rsup|\<ast\>><rsub|1>M<rsub|i
j><rsup|F>n<rsub|1>+[v<rsub|i>]<rsup|\<ast\>><rsub|2>M<rsub|i
j><rsup|F>n<rsub|2><right|]>>>|<row|<cell|>|<cell|=>|<cell|(M<rsup|k>)<rsup|T>\<partial\><rsub|t>\<b-u\>+(S<rsub|\<partial\>x<rsub|1>><rsup|k>)<rsup|T>\<b-v\><rsub|1>+(S<rsub|\<partial\>x<rsub|2>><rsup|k>)<rsup|T>\<b-v\><rsub|2>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><left|[>[v<rsub|F>]<rsup|\<ast\>><rsub|1>(M<rsup|F>)<rsup|T>n<rsub|1>+[v<rsub|F>]<rsup|\<ast\>><rsub|2>(M<rsup|F>)<rsup|T>n<rsub|2><right|]>>>>>
</eqnarray*>
We set <with|mode|math|\<b-v\><rsup|\<ast\>>\<assign\>{\<b-v\>}>. Likewise,
we get
<\eqnarray*>
<tformat|<table|<row|<cell|0>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>v<rsub|t>\<cdot\>\<psi\>-<big|int><rsub|T<rsub|k>>\<nabla\>u\<cdot\>\<psi\>>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>v<rsub|t>\<cdot\>\<psi\>+<big|int><rsub|T<rsub|k>>u\<nabla\>\<cdot\>\<psi\>-<big|int><rsub|T<rsub|k>>\<nabla\>\<cdot\>(u<with|math-font-series|bold|\<psi\>>)>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>v<rsub|t>\<cdot\>\<psi\>+<big|int><rsub|T<rsub|k>>u\<nabla\>\<cdot\>\<psi\>-<big|int><rsub|\<partial\>T<rsub|k>>u\<psi\>\<cdot\>n>>|<row|<cell|>|<cell|\<approx\>>|<cell|<big|int><rsub|T<rsub|k>>v<rsub|t>\<cdot\>\<psi\>+<big|int><rsub|T<rsub|k>>u\<nabla\>\<cdot\>\<psi\>-<big|int><rsub|\<partial\>T<rsub|k>>u<rsup|\<ast\>>\<psi\>\<cdot\>n>>|<row|<cell|>|<cell|\<approx\>>|<cell|<choice|<tformat|<table|<row|<cell|<with|math-display|true|<big|int><rsub|T<rsub|k>><matrix|<tformat|<table|<row|<cell|\<partial\><rsub|t>[v<rsub|i>]<rsub|1>l<rsub|i>>>|<row|<cell|\<partial\><rsub|t>[v<rsub|i>]<rsub|2>l<rsub|i>>>>>>\<cdot\><matrix|<tformat|<table|<row|<cell|l<rsub|j>>>|<row|<cell|0>>>>>+<big|int><rsub|T<rsub|k>>u<rsub|i>l<rsub|i>\<nabla\>\<cdot\><matrix|<tformat|<table|<row|<cell|l<rsub|j>>>|<row|<cell|0>>>>>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><big|int><rsub|F><matrix|<tformat|<table|<row|<cell|u<rsub|i><rsup|\<ast\>>l<rsub|i><rsup|F>l<rsub|j><rsup|F>>>|<row|<cell|u<rsub|i><rsup|\<ast\>>l<rsub|i><rsup|F>0>>>>>\<cdot\>n>>>|<row|<cell|<with|math-display|true|<big|int><rsub|T<rsub|k>><matrix|<tformat|<table|<row|<cell|\<partial\><rsub|t>[v<rsub|i>]<rsub|1>l<rsub|i>>>|<row|<cell|\<partial\><rsub|t>[v<rsub|i>]<rsub|2>l<rsub|i>>>>>>\<cdot\><matrix|<tformat|<table|<row|<cell|0>>|<row|<cell|l<rsub|j>>>>>>+<big|int><rsub|T<rsub|k>>u<rsub|i>l<rsub|i>\<nabla\>\<cdot\><matrix|<tformat|<table|<row|<cell|0>>|<row|<cell|l<rsub|j>>>>>>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><big|int><rsub|F><matrix|<tformat|<table|<row|<cell|u<rsub|i><rsup|\<ast\>>l<rsub|i><rsup|F>0>>|<row|<cell|u<rsub|i><rsup|\<ast\>>l<rsub|i><rsup|F>l<rsub|j><rsup|F>>>>>>\<cdot\>n>>>>>>>>|<row|<cell|>|<cell|=>|<cell|<choice|<tformat|<table|<row|<cell|<with|math-display|true|<big|int><rsub|T<rsub|k>>\<partial\><rsub|t>[v<rsub|i>]<rsub|1>l<rsub|i>l<rsub|j>+<big|int><rsub|T<rsub|k>>u<rsub|i>l<rsub|i>\<partial\><rsub|x<rsub|1>>l<rsub|j>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><big|int><rsub|F>u<rsub|i><rsup|\<ast\>>l<rsub|i><rsup|F>l<rsub|j><rsup|F>n<rsub|x>>>>|<row|<cell|<with|math-display|true|<big|int><rsub|T<rsub|k>>\<partial\><rsub|t>[v<rsub|i>]<rsub|2>l<rsub|i>l<rsub|j>+<big|int><rsub|T<rsub|k>>u<rsub|i>l<rsub|i>\<partial\><rsub|x<rsub|2>>l<rsub|j>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><big|int><rsub|F>u<rsub|i><rsup|\<ast\>>l<rsub|i><rsup|F>l<rsub|j><rsup|F>n<rsub|y>>>>>>>>>|<row|<cell|>|<cell|=>|<cell|<choice|<tformat|<table|<row|<cell|<with|math-display|true|\<partial\><rsub|t>[v<rsub|i>]<rsub|1>M<rsub|i
j><rsup|k>+u<rsub|i>S<rsup|k><rsub|i j,\<partial\>x<rsub|2>>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>>u<rsub|i><rsup|\<ast\>>M<rsub|i
j><rsup|F>n<rsub|x>>>>|<row|<cell|<with|math-display|true|\<partial\><rsub|t>[v<rsub|i>]<rsub|2>M<rsub|i
j><rsup|k>+u<rsub|i>S<rsup|k><rsub|i j,\<partial\>x<rsub|2>>-<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>>u<rsub|i><rsup|\<ast\>>M<rsub|i
j><rsup|F>n<rsub|y>>>>>>>>>>>
</eqnarray*>
and we set <with|mode|math|u<rsup|\<ast\>>\<assign\>{u}>.
<subsubsection|Strong DG>
First equation:
<\eqnarray*>
<tformat|<table|<row|<cell|0>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|D<rsub|k>>v\<cdot\>\<nabla\>\<varphi\>-<big|int><rsub|\<partial\>T<rsub|k>>v<rsup|\<ast\>>\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>(\<nabla\>\<cdot\>v)\<varphi\>+<big|int><rsub|\<partial\>T<rsub|k>>(v-v<rsup|\<ast\>>)\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>(\<nabla\>\<cdot\>v)\<varphi\>+<frac|1|2><big|int><rsub|\<partial\>T<rsub|k>>[v]\<varphi\>\<cdot\>n>>>>
</eqnarray*>
Second equation:
<\eqnarray*>
<tformat|<table|<row|<cell|0>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>v<rsub|t>\<cdot\>\<psi\>+<big|int><rsub|T<rsub|k>>u\<nabla\>\<cdot\>\<psi\>-<big|int><rsub|\<partial\>T<rsub|k>>u<rsup|\<ast\>>\<psi\>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>v<rsub|t>\<cdot\>\<psi\>-<big|int><rsub|T<rsub|k>>(\<nabla\>u)\<psi\>+<big|int><rsub|\<partial\>T<rsub|k>>(u-u<rsup|\<ast\>>)\<psi\>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>v<rsub|t>\<cdot\>\<psi\>-<big|int><rsub|T<rsub|k>>(\<nabla\>u)\<psi\>+<frac|1|2><big|int><rsub|\<partial\>T<rsub|k>>[u]\<psi\>\<cdot\>n>>>>
</eqnarray*>
<subsection|Heat Equation>
Begin with
<\eqnarray*>
<tformat|<table|<row|<cell|u<rsub|t>-\<nabla\>\<cdot\>(<sqrt|a>q)>|<cell|=>|<cell|0,>>|<row|<cell|q-<sqrt|a>\<nabla\>u>|<cell|=>|<cell|0,>>>>
</eqnarray*>
so
<\eqnarray*>
<tformat|<table|<row|<cell|<big|int><rsub|T>[u<rsub|t>-\<nabla\>\<cdot\>(<sqrt|a>q)]\<varphi\>>|<cell|=>|<cell|0>>|<row|<cell|<big|int><rsub|T>u<rsub|t>\<varphi\>+<big|int><rsub|T><sqrt|a>*q\<cdot\>\<nabla\>\<varphi\>-<big|int><rsub|T>\<nabla\>\<cdot\>(<sqrt|a>*q*\<varphi\>)>|<cell|=>|<cell|0>>|<row|<cell|<big|int><rsub|T>u<rsub|t>\<varphi\>+<big|int><rsub|T><sqrt|a>*q\<cdot\>\<nabla\>\<varphi\>-<big|int><rsub|\<partial\>T><sqrt|a>*q*\<varphi\>\<cdot\>n>|<cell|=>|<cell|0>>|<row|<cell|<big|int><rsub|T>u<rsub|t>\<varphi\>+<big|int><rsub|T><sqrt|a>*q\<cdot\>\<nabla\>\<varphi\>-<big|int><rsub|\<partial\>T>(<sqrt|a>*q)<rsup|\<ast\>>\<varphi\>\<cdot\>n>|<cell|=>|<cell|0,>>>>
</eqnarray*>
and
<\eqnarray*>
<tformat|<table|<row|<cell|<big|int><rsub|T>[q-<sqrt|a>\<nabla\>u]\<cdot\>\<psi\>>|<cell|=>|<cell|0>>|<row|<cell|<big|int><rsub|T>q\<cdot\>\<psi\>-<big|int><rsub|T><sqrt|a>\<psi\>\<cdot\>(\<nabla\>u)>|<cell|=>|<cell|0>>|<row|<cell|<big|int><rsub|T>q\<cdot\>\<psi\>+<big|int><rsub|T>\<nabla\>\<cdot\>(<sqrt|a>*u)\<psi\>-<big|int><rsub|T>\<nabla\>\<cdot\>(<sqrt|a>*u\<psi\>)>|<cell|=>|<cell|0>>|<row|<cell|<big|int><rsub|T>q\<cdot\>\<psi\>+<big|int><rsub|T>\<nabla\>\<cdot\>(<sqrt|a>*u)\<psi\>-<big|int><rsub|\<partial\>T>n\<cdot\>(<sqrt|a>*u)<rsup|\<ast\>>\<psi\>>|<cell|=>|<cell|0.>>>> </eqnarray*>
<section|Quadrature Rules>
Golub-Welsch recursion:
<\equation*>
p<rsub|n>=(\<alpha\><rsub|n>x+\<beta\><rsub|n>)p<rsub|n-1>-\<gamma\><rsub|n>p<rsub|n-2>
</equation*>
Hesthaven-Warburton ``recursion'':
<\equation*>
x*p<rsub|n>=a<rsub|n>p<rsub|n-1>+b<rsub|n>p<rsub|n>+a<rsub|n+1>p<rsub|n+1>
</equation*>
Solve for <with|mode|math|a<rsub|n>> and <with|mode|math|b<rsub|n>> from
G-W:
<\eqnarray*>
<tformat|<table|<row|<cell|p<rsub|n+1>>|<cell|=>|<cell|\<alpha\><rsub|n+1>x*p<rsub|n>+\<beta\><rsub|n+1>p<rsub|n>-\<gamma\><rsub|n+1>p<rsub|n-1>>>|<row|<cell|\<alpha\><rsub|n+1>x*p<rsub|n>>|<cell|=>|<cell|p<rsub|n+1>-\<beta\><rsub|n+1>p<rsub|n>+\<gamma\><rsub|n+1>p<rsub|n-1>>>|<row|<cell|x*p<rsub|n>>|<cell|=>|<cell|<frac|1|\<alpha\><rsub|n+1>>p<rsub|n+1>-<frac|\<beta\><rsub|n+1>|\<alpha\><rsub|n+1>>p<rsub|n>+<frac|\<gamma\><rsub|n+1>|\<alpha\><rsub|n+1>>p<rsub|n-1>>>|<row|<cell|>|<cell|=>|<cell|<wide*|<frac|\<gamma\><rsub|n+1>|\<alpha\><rsub|n+1>>|\<wide-underbrace\>><rsub|a<rsub|n>>p<rsub|n-1><wide*|-<frac|\<beta\><rsub|n+1>|\<alpha\><rsub|n+1>>|\<wide-underbrace\>><rsub|b<rsub|n>>p<rsub|n>+<wide*|<frac|1|\<alpha\><rsub|n+1>>|\<wide-underbrace\>><rsub|a<rsub|n+1>>p<rsub|n+1>.>>>>
</eqnarray*>
<section|Upwind Fluxes>
<subsection|DIY Upwind Fluxes/Riemann Solvers>
Start with a linear hyperbolic system:
<\equation*>
\<b-u\><rsub|t>+<big|sum><rsub|i>A<rsub|i>\<partial\><rsub|i>\<b-u\>=0
</equation*>
<\itemize>
<item>Pick a unit normal <math|\<b-n\>> of an interface where the states
<math|\<b-u\><rsup|->> and <math|\<b-u\><rsup|+>> meet.
<item>Compute
<\equation*>
A<rsup|\<pm\>>(\<b-n\>)\<assign\><big|sum><rsub|i>n<rsub|i>A<rsub|i><rsup|\<pm\>>
</equation*>
and diagonalize it to <math|D<rsup|\<pm\>>=V*<rsup|\<pm\>>A<rsup|\<pm\>>(\<b-n\>)*(V<rsup|\<pm\>>)<rsup|T>>,
resulting in eigenvalues
<\equation*>
\<lambda\><rsub|1><rsup|\<pm\>>\<leqslant\>\<cdots\>\<leqslant\>\<lambda\><rsub|k><rsup|\<pm\>>\<leqslant\>0\<leqslant\>\<lambda\><rsub|k+1><rsup|\<pm\>>\<leqslant\>\<cdots\>\<leqslant\>\<lambda\><rsub|n><rsup|\<pm\>>.
</equation*>
Call <math|\<b-s\><rsup|\<pm\>>\<assign\>V<rsup|\<pm\>>\<b-u\><rsup|\<pm\>>>.
<item>Consider the space-time diagram of Figure <reference|fig:fluxfan>.
Each eigenvalue is imagined to be the speed of a <em|shock> emanating
from the point under consideration, and each shock separating one
<em|state> from another. The left- and rightmost state are
<math|\<b-s\><rsup|->> and <math|\<b-s\><rsup|+>>, of course. The next
state from the left is labeled <math|\<b-s\><rsup|\<ast\>>>, the one
after that <math|\<b-s\><rsup|\<ast\>\<ast\>>>, up to
<math|\<b-s\><rsup|\<ast\>(n-1)>>. These states are a-priori unknown. In
addition, the fluxes <math|(D\<b-s\>)<rsup|\<ast\>>,\<ldots\>,(D\<b-s\>)<rsup|\<ast\>(n-1)>>
are also unknown: For the moment, there are <math|2(n-1)> unknowns.
<item>The Rankine-Hugoniot condition
<\equation*>
\<lambda\><rsub|i><rsup|\<pm\>>=<frac|(D\<b-s\>)<rsup|*\<ast\>(i)>-(D\<b-s\>)<rsup|*\<ast\>(i-1)>|\<b-s\><rsup|*\<ast\>(i)>-\<b-s\><rsup|*\<ast\>(i-1)>>
</equation*>
has to hold across each such interface. For convenience, set <math|\<b-s\><rsup|\<ast\>(0)>=\<b-s\><rsup|->> and
<math|\<b-s\><rsup|\<ast\>(n)>=\<b-s\><rsup|+>>. This will provide
<math|n> equations.
Naturally, eigenvalues larger than zero (``right-travelling'') use their
<math|\<cdot\><rsup|->> values, whereas left-travelling values smaller
than zero use their <math|\<cdot\><rsup|+>> values.
<item>We label <math|\<lambda\><rsub|k>> that eigenvalue whose state
straddles zero the line of zero speed--that includes the case of
<math|\<lambda\><rsub|k>=0>. If there are several zero eigenvalues (can
that even happen?), pick <math|k> as large as possible. It is assumed
that the number <math|k> is the same on both sides of the interface, and
therefore it carries no <math|\<pm\>> superscript.
<item>The numerical flux, i.e. the solution that we seek, is the flux
<math|(D\<b-s\>)<rsup|\<ast\>(k)>> that straddles zero.
<item>For non-zero-straddling states, the unknowns
<math|(D\<b-s\>)<rsup|\<ast\>(i)>> are found as follows:
<\itemize>
<item>If <math|i\<less\>k>, <math|(D\<b-s\>)<rsup|\<ast\>(i)>=D<rsup|->\<b-s\><rsup|\<ast\>(i)>>.
<item>If <math|i\<gtr\>k>, <math|(D\<b-s\>)<rsup|\<ast\>(i)>=D<rsup|+>\<b-s\><rsup|\<ast\>(i)>>.
</itemize>
This takes <math|n-2> unknowns out of the picture, leaving us with
<math|2(n-1)-(n-2)=n> unknowns.
<item>The case of a zero eigenvalue deserves special mention, even though
its handling is technically no different than before: The states
<math|\<b-s\><rsup|\<ast\>(k-1)>> and <math|\<b-s\><rsup|\<ast\>(k)>>
fall cleanly to the left and right of zero, so we can use
<\eqnarray*>
<tformat|<table|<row|<cell|(D\<b-s\>)<rsup|\<ast\>(k-1)>>|<cell|=>|<cell|D<rsup|->\<b-s\><rsup|\<ast\>(k-1)>,>>|<row|<cell|(D\<b-s\>)<rsup|\<ast\>(k)>>|<cell|=>|<cell|D<rsup|+>\<b-s\><rsup|\<ast\>(k)>.>>>>
</eqnarray*>
Observe that we end up with only <math|n-1> unknowns (namely all the
<math|{\<b-s\><rsup|\<ast\>(i)>}<rsub|i=1><rsup|n-1>>), but we still have
<math|n> equations. It seems, however, that the resulting system of
equations has rank <math|n-1>, and so a unique solution can be
determined. (<with|color|red|Why is that, rigorously?>)
Lastly, observe that the <math|k>th Rankine-Hugoniot condition dictates
<\equation*>
(D\<b-s\>)<rsup|\<ast\>(k-1)>=(D\<b-s\>)<rsup|\<ast\>(k)>,
</equation*>
ensuring consistency across the zero-speed shock.
</itemize>
I believe that the above can also be applied in the case of degenerate
eigenvalues, but haven't checked in detail.
<big-figure|<postscript|fluxfan.fig|9cm|||||>|<label|fig:fluxfan>Space-Time
Diagram of a flux fan.>
<subsection|Comments by Akil>
<subsubsection|Zero eigenvalues>
I'm not entirely sure that the case of zero eigenvalues actually comes out
well. Let's consider the following: adopting your notation above, let's
assume we have <math|n> eigenvalues for each element,
<math|\<lambda\><rsub|i><rsup|\<pm\>>>, <math|i = 1,2,\<ldots\>n>. Let
<math|N>, <math|M>, and <math|P> represent the number of negative, zero,
and positive eigenvalues, respectively. (Recall that we're assuming that <math|sign(\<lambda\><rsub|i><rsup|+>) = sign(\<lambda\><rsub|i><rsup|->)>
for all <math|i>.) We have <math|N+M+P=n>, and we assume at least two of
these numbers are nonzero. (If only one of them is nonzero, upwinding is
trivial.)
\;
We start with the <math|2*(n-1)> unknowns
<math|<left|{><with|math-font-series|bold|s><rsup|\<asterisk\>(i)>,
<left|(>D*<with|math-font-series|bold|s><right|)><rsup|\<asterisk\>(i)><right|}><rsub|i=1><rsup|n-1>>.\
\;
I'm adopting the convention that positive eigenvalues correspond to
right-traveling waves. We have the following Rankine-Hugoniot conditions:
<\equation>
<tabular|<tformat|<cwith|3|3|3|3|cell-halign|l>|<table|<row|<cell|\<lambda\><rsup|-><rsub|i>>|<cell|=>|<cell|<frac|(D\<b-s\>)<rsup|*\<ast\>(i)>-(D\<b-s\>)<rsup|*\<ast\>(i-1)>|\<b-s\><rsup|*\<ast\>(i)>-\<b-s\><rsup|*\<ast\>(i-1)>>,>|<cell|>|<cell|i
= 1, 2, \<ldots\>,N>>|<row|<cell|>|<cell|>|<cell|>|<cell|>|<cell|>>|<row|<cell|0>|<cell|=>|<cell|(D\<b-s\>)<rsup|*\<ast\>(i)>-(D\<b-s\>)<rsup|*\<ast\>(i-1)>,>|<cell|>|<cell|i
= N+1, \<ldots\>,N+Z>>|<row|<cell|>|<cell|>|<cell|>|<cell|>|<cell|>>|<row|<cell|\<lambda\><rsup|+><rsub|i>>|<cell|=>|<cell|<frac|(D\<b-s\>)<rsup|*\<ast\>(i)>-(D\<b-s\>)<rsup|*\<ast\>(i-1)>|\<b-s\><rsup|*\<ast\>(i)>-\<b-s\><rsup|*\<ast\>(i-1)>>,>|<cell|>|<cell|i
= N+Z+1, \<ldots\>,n>>>>><right|}> \ n \ (linear)
equations<label|eq:rh-all>
</equation>
\;
\;
If I understand what we're doing, we consider the unknowns <math|<left|{>
<left|(>D*<with|math-font-series|bold|s><right|)><rsup|\<asterisk\>(i)><right|}><rsub|i=1,
2, \<ldots\><with|math-font-series|bold|N-1>,
<with|math-font-series|bold|N+Z+1>, \<ldots\>, n>> as anisotropic
(directionally-biased) unknowns with explicit connections to the
<math|<with|math-font-series|bold|s><rsup|\<asterisk\>(i)>> given by the
auxiliary equations
<\equation>
<tabular|<tformat|<table|<row|<cell|(D\<b-s\>)<rsup|*\<ast\>(i)> =
D<rsup|->\<b-s\><rsup|*\<ast\>(i)>>|<cell|>|<cell|i = 1, 2,
\<ldots\>,N-1>>|<row|<cell|>|<cell|>|<cell|>>|<row|<cell|(D\<b-s\>)<rsup|*\<ast\>(i)>
= D<rsup|+>\<b-s\><rsup|*\<ast\>(i)>>|<cell|>|<cell|i = N+Z+1, \<ldots\>,
n>>>>><right|}> \ N+P-2 \ \ (linear) equations<label|eq:rh-aux>
</equation>
\;
However, I don't see what we do for the unknowns which make contact with
the zero eigenvalue(s). You mention that if <math|Z=1> with
<math|\<lambda\><rsub|k>=0> we would have\
<\equation*>
<tabular|<tformat|<table|<row|<cell|(D\<b-s\>)<rsup|\<ast\>(k)>>|<cell|<above|=|<with|mode|text|NO!>>>|<cell|D<rsup|+>\<b-s\><rsup|\<ast\>(k)>.>>>>>
</equation*>
But then I'd also claim that\
<\equation*>
<tabular|<tformat|<table|<row|<cell|(D\<b-s\>)<rsup|\<ast\>(k-1)>>|<cell|<above|=|<with|mode|text|NO!>>>|<cell|D<rsup|->\<b-s\><rsup|\<ast\>(k-1)>.>>>>>
</equation*>
I cannot think of a reason why we'd treat
<with|mode|math|(D\<b-s\>)<rsup|\<ast\>(k-1)>> any differently from
<with|mode|math|(D\<b-s\>)<rsup|\<ast\>(k)>> since I don't think it's the
case that <with|mode|math|(D\<b-s\>)<rsup|\<ast\>(k)>> `belongs' to
<math|\<lambda\><rsub|k>=0> in any way. If <math|Z=0>, then indeed
(<reference|eq:rh-aux>) gives us <math|N+P-2 = n-2> equations to add to the
<math|n> equations in (<reference|eq:rh-all>) giving <math|2*n-2> equations
to couple the <math|2(n-1)> unknowns. But if <math|Z=1>, then
(<reference|eq:rh-aux>) gives us <math|N+P-2=n-3> equations to couple with
the <math|n> from (<reference|eq:rh-all>). That's <math|2*n-3> equations for <math|2*n-2> unknowns. Conversely, if instead we adopted the convention
<\equation*>
<tabular|<tformat|<table|<row|<cell|(D\<b-s\>)<rsup|\<ast\>(k)>>|<cell|<above|=|<with|mode|text|YES!>>>|<cell|D<rsup|+>\<b-s\><rsup|\<ast\>(k)>>>|<row|<cell|>|<cell|>|<cell|>>|<row|<cell|(D\<b-s\>)<rsup|\<ast\>(k-1)>>|<cell|<above|=|<with|mode|text|YES!>>>|<cell|D<rsup|+>\<b-s\><rsup|\<ast\>(k)>,>>>>>
</equation*>
then we'd have 2 equations plus <math|N+P-2=n-3> from
(<reference|eq:rh-aux>) plus <math|n> from (<reference|eq:rh-all>) giving
us <math|2*n-1> equations for <math|2*n-2> unknowns.\
\;
In general, if we don't adopt any directional bias for states adjacent to
the zero eigenvalue, then for <math|Z\<gtr\>0> zero eigenvalues, we'd have
<math|2*n-2-Z> equations for <math|2*n-2> unknowns, and if we do adopt
directional bias for states adjacent to the zero eigenvalue, then we'd have
<math|2n-Z> equations for <math|2n-2> unknowns (If
<math|\<lambda\><rsub|k>=\<lambda\><rsub|k+1>=0>, I don't see how you can
adopt a directional bias for <with|mode|math|(D\<b-s\>)<rsup|\<ast\>(k)>>.)
\;
You make the comment above that ``We simply regard
<math|(D\<b-s\>)<rsup|\<ast\>(k)>> as another unknown of the system, which,
in addition to <math|{\<b-s\><rsup|*\<ast\>(i)>}<rsub|1><rsup|n-1>>, brings
us to <math|n> unknowns.''
\;
If I understand this correctly, this indeed does make <math|n> equations
and <math|n> unknowns, but you treat <with|mode|math|(D\<b-s\>)<rsup|\<ast\>(k)>>
differently than <with|mode|math|(D\<b-s\>)<rsup|\<ast\>(k+1)>> (non
directionally-biased vs. directionally biased) and, again, I don't see how
you can consistently make this choice given that, to me,
<with|mode|math|(D\<b-s\>)<rsup|\<ast\>(k+1)>> is no more directionally
unbiased than <with|mode|math|(D\<b-s\>)<rsup|\<ast\>(k)>>.\
<subsubsection|Degenerate eigenvalues>
\;
Another comment: degenerate eigenvalues should be ok: let's suppose that
<math|\<lambda\><rsub|k>=\<lambda\><rsub|k+1>>. Then there is an
intermediate state <math|<with|math-font-series|bold|s><rsup|\<asterisk\>(k)>>
separating <math|<with|math-font-series|bold|s><rsup|\<asterisk\>(k-1)>>
and <math|<with|math-font-series|bold|s><rsup|\<asterisk\>(k+1)>>. This
will give us\
<\equation*>
<tabular|<tformat|<table|<row|<cell|\<lambda\><rsub|k>=<frac|(D\<b-s\>)<rsup|*\<ast\>(k+1)>-(D\<b-s\>)<rsup|*\<ast\>(k)>|\<b-s\><rsup|*\<ast\>(k+1)>-\<b-s\><rsup|*\<ast\>(k)>>>>|<row|<cell|>>|<row|<cell|\<lambda\><rsub|k>=<frac|(D\<b-s\>)<rsup|*\<ast\>(k)>-(D\<b-s\>)<rsup|*\<ast\>(k-1)>|\<b-s\><rsup|*\<ast\>(k)>-\<b-s\><rsup|*\<ast\>(k-1)>>>>>>>
</equation*>
Using these equations to eliminate the unknowns
<with|mode|math|(D\<b-s\>)<rsup|*\<ast\>(k)>> and
<with|mode|math|\<b-s\><rsup|*\<ast\>(k)>> gives us the Rankine-Hugoniot
condition relating <with|mode|math|\<b-s\><rsup|*\<ast\>(k-1)>> to
<with|mode|math|\<b-s\><rsup|*\<ast\>(k+1)>>:
<\equation*>
\<lambda\><rsub|k>=<frac|(D\<b-s\>)<rsup|*\<ast\>(k+1)>-(D\<b-s\>)<rsup|*\<ast\>(k-1)>|\<b-s\><rsup|*\<ast\>(k+1)>-\<b-s\><rsup|*\<ast\>(k-1)>>.
</equation*>
This is the same condition one would get if we pretended that states
<with|mode|math|(D\<b-s\>)<rsup|*\<ast\>(k)>> and
<with|mode|math|\<b-s\><rsup|*\<ast\>(k)>>, the states straddled by the
shocks of the same speed, did not exist and merrily went about our
upwinding.\
\;
Note that this also means that a <math|Z>-fold degeneracy in the zero
eigenvalue isn't a problem either, as long as we can sort out the simple
zero eigenvalue case.
<subsection|More random comments by Akil>
We seem to be pretty confident that if there does not exist any 0
eigenvalue, then everything is fine. Kittens meow, rainbows shine, and the
world is happy. But we've also established that we don't know what's going
on when <math|\<exists\> \ \ \<lambda\>=0>. However, let's boil things down
to the simplest case we can. We know degenerate eigenvalues pose no problem
(see my comment last time), and multiple simple (+) or (<math|->)
eigenvalues are likewise fine. So let's consider this model case:
<\equation*>
A<rsup|\<pm\>>(<wide|<with|math-font-series|bold|n>|^>):<with|math-font|Bbb*|R><rsup|3>\<rightarrow\><with|math-font|Bbb*|R><rsup|3>,
\ \ \<lambda\><rsup|\<pm\>><rsub|1>\<less\>0,
\<lambda\><rsup|\<pm\>><rsub|2>=0, \ \ \ and
\ \ \ \<lambda\><rsup|\<pm\>><rsub|3>\<gtr\>0.
</equation*>
We assume the greatest generality possible here:
<math|\<lambda\><rsub|\<cdot\>><rsup|\<pm\>>> can be different magnitudes,
but they have the same sign. In line with a strong hyperbolic system, we
assume that <math|A<rsup|\<pm\>>> is diagonalizable as
<\equation*>
A<rsup|\<pm\>>=V<rsup|\<pm\>>\<Lambda\><rsup|\<pm\>><left|(>V<rsup|\<pm\>><right|)><rsup|T>,
</equation*>
where <math|\<Lambda\><rsup|\<pm\>>> are diagonal matrices containing the
<math|\<lambda\><rsup|\<pm\>><rsub|i>>. We introduce the eigenvectors
<\equation*>
(V<rsup|\<pm\>>) =<left|(><tabular|<tformat|<table|<row|<cell|<with|math-font-series|bold|v><rsub|1><rsup|\<pm\>>>|<cell|<with|math-font-series|bold|v><rsub|2><rsup|\<pm\>>>|<cell|<with|math-font-series|bold|v><rsub|3><rsup|\<pm\>>>>>>><right|)>
</equation*>
Here, all eigenvalues are simple, and <math|<left|{><with|math-font-series|bold|v><rsub|i><rsup|\<pm\>><right|}><rsub|i=1><rsup|3>>
span <math|<with|math-font|Bbb*|R><rsup|3>> for each (+) or (<math|->).
<with|font-series|bold|Note:> there is no reason to assume that, e.g.
<math|<left|{><with|math-font-series|bold|v><rsub|1><rsup|+>,
<with|math-font-series|bold|v><rsub|2><rsup|->,
<with|math-font-series|bold|v><rsub|3><rsup|-><right|}>> span
<math|<with|math-font|Bbb*|R><rsup|3>>. Or we shall not assume so here, at
least. Although, to be fair, if this is not the case, then we can have the
case of a genuine shock forming, which is kind of ridiculous. Of course, we
have that <math|\<lambda\><rsub|i><rsup|\<pm\>>> is associated with
<math|<with|math-font-series|bold|v><rsub|i><rsup|\<pm\>>>. We also assume
that\
<\equation*>
<tabular|<tformat|<table|<row|<cell|<with|math-font-series|bold|u><rsup|->=<big|sum><rsub|i=1><rsup|3>\<mu\><rsup|-><rsub|i><with|math-font-series|bold|v><rsub|i><rsup|->>|<cell|>|<cell|<with|math-font-series|bold|u><rsup|+>=<big|sum><rsub|i=1><rsup|3>\<mu\><rsub|i><rsup|+><with|math-font-series|bold|v><rsub|i><rsup|+>>>>>>,
</equation*>
which we can do since the <math|<with|math-font-series|bold|v><rsub|i>>'s
span. And of course, I just realized that the
<math|\<mu\><rsub|i><rsup|\<pm\>>> are the entries in
<math|<with|math-font-series|bold|s><rsup|\<pm\>><rsub|i>>. Go math.\
\;
We introduce the intermediate states <math|<with|math-font-series|bold|u><rsup|\<asterisk\>>>
and <math|<with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>>,
with their appropriate fluxes <math|(A<with|math-font-series|bold|u>)<rsup|\<asterisk\>>>
and <math|(A<with|math-font-series|bold|u>)<rsup|\<asterisk\>\<asterisk\>>>.
All four of these, at the moment, are unknown.\
\;
We note that at the end of the day, we want some flux\
<\equation*>
(A<with|math-font-series|bold|u>)<rsup|#>
</equation*>
(yes that's horrible notation; leave me alone! wahhh!)
\;
All of our algorithms will attempt to find this quantity. Now, we consider
three possible algorithms for upwind flux computations.\
<subsubsection|Stoopid upwinding>
The motivation here is simple: the 255 view of upwinding is the following:
the evolution equation for <math|u<rsup|->> should not allow boundary
conditions that impose <math|<with|math-font-series|bold|v><rsub|1><rsup|->>;
similarly, the evolution equation for <math|u<rsup|+>> should not allow
boundary conditions that impose <math|<with|math-font-series|bold|v><rsub|3><rsup|+>>.
This is algorithmically consistent with vanilla upwinding
<math|u<rsub|t>=a*u<rsub|x>>. The algorithm then is the following:
<\enumerate>
<item>We simply set <math|<with|math-font-series|bold|u><rsup|\<asterisk\>>=<with|math-font-series|bold|u><rsup|->-\<mu\><rsup|-><rsub|2>*<with|math-font-series|bold|v><rsup|-><rsub|2>-\<mu\><rsup|-><rsub|1><with|math-font-series|bold|v><rsup|-><rsub|1>=\<mu\><rsup|-><rsub|3><with|math-font-series|bold|v><rsup|-><rsub|3>>.
Also set <math|<with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>=<with|math-font-series|bold|u><rsup|+>-\<mu\><rsup|+><rsub|2>*<with|math-font-series|bold|v><rsup|+><rsub|2>-\<mu\><rsup|+><rsub|3><with|math-font-series|bold|v><rsup|+><rsub|3>=\<mu\><rsub|1><rsup|+><with|math-font-series|bold|v><rsup|+><rsub|1>>.\
<item>Set <math|(A<with|math-font-series|bold|u>)<rsup|\<asterisk\>>=A<rsup|-><with|math-font-series|bold|u><rsup|\<asterisk\>>>
and <math|(A*<with|math-font-series|bold|u>)<rsup|\<asterisk\>\<asterisk\>>=A<rsup|+>*<with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>>.\
<item>Use the flux <math|(A<with|math-font-series|bold|s>)<rsup|#>=(A<with|math-font-series|bold|u>)<rsup|\<asterisk\>>+(A*<with|math-font-series|bold|u>)<rsup|\<asterisk\>\<asterisk\>>=\<mu\><rsup|-><rsub|3><with|math-font-series|bold|v><rsup|-><rsub|3>+\<mu\><rsub|1><rsup|+><with|math-font-series|bold|v><rsup|+><rsub|1>>.
</enumerate>
\;
Advantages of this flux:
<\itemize>
<item>It's easy
<item>It's very intuitive
<item>It's consistent with usual upwinding if <math|A<rsup|+>=A<rsup|->>.
</itemize>
Disadvantages:
<\itemize>
<item>There is no motivation for it other than, ``it's upwinding''.
<item>Where'r the R-H conditions? Actually, these are
<with|font-series|bold|almost> exactly those conditions: see below....
</itemize>
<subsubsection|R-H Conditions with Directional Bias>
In this formulation, we impose the Rankine-Hugoniot conditions, which in
this case read
<\equation*>
<tabular|<tformat|<cwith|3|3|1|1|cell-halign|r>|<table|<row|<cell|\<lambda\><rsub|1><rsup|->(<with|math-font-series|bold|u><rsup|\<asterisk\>>-<with|math-font-series|bold|u><rsup|->)>|<cell|=>|<cell|<left|(>A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>>-A<rsup|-><with|math-font-series|bold|u><rsup|->>>|<row|<cell|>|<cell|>|<cell|>>|<row|<cell|(A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>>>|<cell|=>|<cell|(A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>\<asterisk\>>>>|<row|<cell|>|<cell|>|<cell|>>|<row|<cell|\<lambda\><rsub|3><rsup|+>(<with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>-<with|math-font-series|bold|u><rsup|+>)>|<cell|=>|<cell|<left|(>A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>\<asterisk\>>-A<rsup|+><with|math-font-series|bold|u><rsup|+>>>>>>
</equation*>
\;
However, we also impose `directional bias', which means that intermediate
states which are adjacent to, but not overlapping, the zero-speed shock
inherit the operators from their respective sides. In math, this reads
<\equation*>
<tabular|<tformat|<table|<row|<cell|(A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>>=A<rsup|-><with|math-font-series|bold|u><rsup|\<asterisk\>>>|<cell|>|<cell|(A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>\<asterisk\>>=A<rsup|+><with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>.>>>>> </equation*>
This is very interesting set of conditions: it seemingly over-determines
the system, and the R-H conditions become\
<\equation*>
<tabular|<tformat|<cwith|3|3|1|1|cell-halign|r>|<table|<row|<cell|<left|(>A<rsup|->-\<lambda\><rsub|1><rsup|->I<right|)><left|(><with|math-font-series|bold|u><rsup|\<asterisk\>>-<with|math-font-series|bold|u><rsup|-><right|)>>|<cell|=>|<cell|0>>|<row|<cell|>|<cell|>|<cell|>>|<row|<cell|A<rsup|-><with|math-font-series|bold|u><rsup|\<asterisk\>>>|<cell|=>|<cell|A<rsup|+><with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>>>|<row|<cell|>|<cell|>|<cell|>>|<row|<cell|<left|(>A<rsup|+>-\<lambda\><rsub|3><rsup|+>I<right|)><left|(><with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>-<with|math-font-series|bold|u><rsup|+><right|)>>|<cell|=>|<cell|0>>>>>
</equation*>
\;
This means that the quantity <with|mode|math|<with|math-font-series|bold|u><rsup|\<asterisk\>>-<with|math-font-series|bold|u><rsup|->>
*must* be a multiple of <math|<with|math-font-series|bold|v><rsup|-><rsub|1>>.
I.e.,\
<\equation*>
<with|math-font-series|bold|u><rsup|\<asterisk\>>=\<nu\><rsub|1><with|math-font-series|bold|v><rsub|1><rsup|->+\<mu\><rsub|2><rsup|-><with|math-font-series|bold|v><rsub|2><rsup|->+\<mu\><rsub|3><rsup|-><with|math-font-series|bold|v><rsub|3><rsup|->,
</equation*>
for some unknown scalar <math|\<nu\><rsub|1>>. A similar argument leads us
to\
<\equation*>
<with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>=\<mu\><rsup|+><rsub|1><with|math-font-series|bold|v><rsub|1><rsup|+>+\<mu\><rsub|2><rsup|+><with|math-font-series|bold|v><rsub|2><rsup|+>+\<nu\><rsub|3><with|math-font-series|bold|v><rsub|3><rsup|+>.
</equation*>
It is worth stressing that our unknown vectors
<math|<with|math-font-series|bold|u><rsup|\<asterisk\>>> and
<math|<with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>>,
previously a full 6 unknowns, are now merely two unknowns,
<math|\<nu\><rsub|1>> and <math|\<nu\><rsub|3>>. This will be true for any
size of the matrices <math|A<rsup|\<pm\>>>.\
\;
<with|font-shape|italic|Note that this has a strikingly similar flavor to
`stoopid upwinding': in that case, we would have <math|\<nu\><rsub|1>=0>
and <math|\<nu\><rsub|3>=0>. If we impose directional bias, `Stoopid
upwinding' automatically satisfies all but one of the R-H conditions: the
one at <math|\<lambda\>=0>.>
\;
Of course, we are still overdetermined: we have one R-H condition left:
three equations, two unknowns. What's left? The <math|\<lambda\>=0>
conditions boils down to\
<\equation*>
\<nu\><rsub|1>\<lambda\><rsub|1><rsup|->*<with|math-font-series|bold|v><rsub|1><rsup|->-\<nu\><rsub|3>*\<lambda\><rsub|3><rsup|+><with|math-font-series|bold|v><rsub|3><rsup|+>=\<mu\><rsub|1><rsup|+>\<lambda\><rsub|1><rsup|+><with|math-font-series|bold|v><rsub|1><rsup|+>-\<mu\><rsub|3><rsup|->\<lambda\><rsub|3><rsup|-><with|math-font-series|bold|v><rsub|3><rsup|->,
</equation*>
where the unknowns are <math|\<nu\><rsub|1>> and <math|\<nu\><rsub|3>>.
Note that we *need* this equation to be satisfied independent of the
unknowns <math|\<mu\><rsub|1><rsup|+>> and <math|\<mu\><rsub|3><rsup|->>,
which are dependent on the incoming waves
<math|<with|math-font-series|bold|u><rsup|\<pm\>>>, over which we have no
control. This means that we require that
<math|span<left|{><with|math-font-series|bold|v><rsub|1><rsup|+>,
<with|math-font-series|bold|v><rsub|3><rsup|-><right|}>\<subset\>span<left|{><with|math-font-series|bold|v><rsub|1><rsup|->,<with|math-font-series|bold|v><rsub|3><rsup|+><right|}>>.
Note that since <math|<with|math-font-series|bold|v><rsub|1><rsup|+>\<perp\><with|math-font-series|bold|v><rsub|3><rsup|+>>,
and <math|<with|math-font-series|bold|v><rsub|1><rsup|->\<perp\><with|math-font-series|bold|v><rsub|3><rsup|->>,
then we must have that <math|<with|math-font-series|bold|v><rsub|1><rsup|+>\<parallel\><with|math-font-series|bold|v><rsub|1><rsup|+>>
and <math|<with|math-font-series|bold|v><rsub|3><rsup|->\<parallel\><with|math-font-series|bold|v><rsub|3><rsup|->>.
Then (and only then!) will we be able to solve for the fluxes we wish for
any values of <math|\<mu\><rsup|\<pm\>><rsub|\<cdot\>>>.\
\;
In this simple case, our condition states that
<math|<with|math-font-series|bold|v><rsub|1><rsup|+>=c<rsub|1>*<with|math-font-series|bold|v><rsub|1><rsup|->> and <math|<with|math-font-series|bold|v><rsub|3><rsup|+>=c<rsub|3>*<with|math-font-series|bold|v><rsub|3><rsup|->>.
I.e., that left traveling waves stay left-traveling waves and
right-traveling waves do the same. Then we have that\
<\equation*>
<tabular|<tformat|<table|<row|<cell|\<nu\><rsub|1>=c<rsub|1>*\<mu\><rsub|1><rsup|+>*<frac|\<lambda\><rsub|1><rsup|+>|\<lambda\><rsub|1><rsup|->>>|<cell|>|<cell|\<nu\><rsub|3>=\<mu\><rsub|3><rsup|-><frac|\<lambda\><rsub|3><rsup|->|c<rsub|3>*\<lambda\><rsub|3><rsup|+>>>>>>>,
</equation*>
which uniquely determines the flux as\
<\equation*>
<left|(>A<with|math-font-series|bold|u><right|)><rsup|#>=A<rsup|-><with|math-font-series|bold|u><rsup|\<asterisk\>>=A<rsup|+><with|math-font-series|bold|u><rsup|\<asterisk\>\<asterisk\>>=\<nu\><rsub|1>\<lambda\><rsub|1><rsup|-><with|math-font-series|bold|v><rsub|1><rsup|->+\<mu\><rsub|3><rsup|->\<lambda\><rsub|3><rsup|-><with|math-font-series|bold|v><rsub|3><rsup|->=\<mu\><rsup|+><rsub|1>\<lambda\><rsub|1><rsup|+><with|math-font-series|bold|v><rsub|1><rsup|+>+\<nu\><rsub|3>\<lambda\><rsub|3><rsup|+><with|math-font-series|bold|v><rsub|3><rsup|+>.
</equation*>
\;
The above considerations lead us to the following theorem:
<\theorem>
(Well-posedness of Rankine-Hugoniot directionally-biased upwinding)
Let <math|A<rsup|\<pm\>>(<wide|<with|math-font-series|bold|n>|^>)>
represent the (discontinuous) operator of a multidimensional strongly
hyperbolic linear wave equation at an element interface. Assume
diagonalizations of the form\
<\equation*>
<tabular|<tformat|<table|<row|<cell|A<rsup|\<pm\>>=V<rsup|\<pm\>>\<Lambda\><rsup|\<pm\>><left|(>V<rsup|\<pm\>><right|)><rsup|T>,>>|<row|<cell|>>|<row|<cell|V<rsup|\<pm\>>=<left|(><tabular|<tformat|<table|<row|<cell|<with|math-font-series|bold|v><rsup|\<pm\>><rsub|1>>|<cell|<with|math-font-series|bold|v><rsub|2><rsup|\<pm\>>>|<cell|\<cdots\>>|<cell|<with|math-font-series|bold|v><rsub|n><rsup|\<pm\>>>>>>><right|)>>>>>>
</equation*>
\;
Let the evolving <math|<with|math-font-series|bold|u>\<in\><with|math-font|Bbb*|R><rsup|n>>.
The known values of <math|<with|math-font-series|bold|u>> at the element
interfaces admit the decomposition\
<\equation*>
<tabular|<tformat|<table|<row|<cell|<with|math-font-series|bold|u><rsup|->=<big|sum><rsub|i=1><rsup|n>\<mu\><rsub|i><rsup|-><with|math-font-series|bold|v><rsub|i><rsup|->>|<cell|>|<cell|<with|math-font-series|bold|u><rsup|+>=<big|sum><rsub|i=1><rsup|n>\<mu\><rsub|i><rsup|+><with|math-font-series|bold|v><rsub|i><rsup|+>.>>>>>
</equation*>
Assume that each diagonal matrix <math|\<Lambda\><rsup|\<pm\>>> contains
eigenvalues <math|\<lambda\><rsub|i><rsup|\<pm\>>> such that the
following quantities are well-defined:
<\equation*>
<tabular|<tformat|<table|<row|<cell|<with|math-font|cal|P><rsup|+>=<left|{>i:
\<lambda\><rsup|+><rsub|i>\<gtr\>0<right|}>>|<cell|>|<cell|<with|math-font|cal|N><rsup|+>=<left|{>i:
\<lambda\><rsup|+><rsub|i>\<less\>0<right|}>>|<cell|>|<cell|<with|math-font|cal|Z><rsup|+>=<left|{>i:
\<lambda\><rsup|+><rsub|i>=0<right|}>>>|<row|<cell|>|<cell|>|<cell|>|<cell|>|<cell|>>|<row|<cell|<with|math-font|cal|P><rsup|->=<left|{>i:
\<lambda\><rsup|-><rsub|i>\<gtr\>0<right|}>>|<cell|>|<cell|<with|math-font|cal|N><rsup|->=<left|{>i:
\<lambda\><rsup|-><rsub|i>\<less\>0<right|}>>|<cell|>|<cell|<with|math-font|cal|Z><rsup|->=<left|{>i:
\<lambda\><rsup|-><rsub|i>=0<right|}>>>|<row|<cell|>|<cell|>|<cell|>|<cell|>|<cell|>>|<row|<cell|P
=<left|\|><with|math-font|cal|P><rsup|+><right|\|>=<left|\|><with|math-font|cal|P><rsup|-><right|\|>>|<cell|>|<cell|N
=<left|\|><with|math-font|cal|N><rsup|+><right|\|>=<left|\|><with|math-font|cal|N><rsup|-><right|\|>>|<cell|>|<cell|Z
=<left|\|><with|math-font|cal|Z><rsup|+><right|\|>=<left|\|><with|math-font|cal|Z><rsup|-><right|\|>.>>>>>
</equation*>
Note that this gives us <math|P+Z+N=n>. Associate with each eigenvalue
<math|\<lambda\><rsup|\<pm\>><rsub|i>> the eigenvector
<math|<with|math-font-series|bold|v><rsub|i><rsup|\<pm\>>>. Define the
following vector subspaces of <math|<with|math-font|Bbb*|R><rsup|n>>:
<\equation*>
<tabular|<tformat|<table|<row|<cell|<with|math-font|cal|W><rsub|<with|math-font|cal|P>><rsup|+>=span<left|{><with|math-font-series|bold|v><rsub|i><rsup|+>:i\<in\><with|math-font|cal|P><rsup|+><right|}>>|<cell|>|<cell|<with|math-font|cal|W><rsub|<with|math-font|cal|N>><rsup|+>=span<left|{><with|math-font-series|bold|v><rsub|i><rsup|+>:i\<in\><with|math-font|cal|N><rsup|+><right|}>>|<cell|>|<cell|<with|math-font|cal|W><rsub|<with|math-font|cal|Z>><rsup|+>=span<left|{><with|math-font-series|bold|v><rsub|i><rsup|+>:i\<in\><with|math-font|cal|Z><rsup|+><right|}>>>|<row|<cell|>|<cell|>|<cell|>|<cell|>|<cell|>>|<row|<cell|<with|math-font|cal|W><rsub|<with|math-font|cal|P>><rsup|->=span<left|{><with|math-font-series|bold|v><rsub|i><rsup|->:i\<in\><with|math-font|cal|P><rsup|-><right|}>>|<cell|>|<cell|<with|math-font|cal|W><rsub|<with|math-font|cal|N>><rsup|->=span<left|{><with|math-font-series|bold|v><rsub|i><rsup|->:i\<in\><with|math-font|cal|N><rsup|-><right|}>>|<cell|>|<cell|<with|math-font|cal|W><rsub|<with|math-font|cal|Z>><rsup|->=span<left|{><with|math-font-series|bold|v><rsub|i><rsup|->:i\<in\><with|math-font|cal|Z><rsup|-><right|}>>>>>>
</equation*>
\;
Define the <math|(n-1)> intermediate states
<math|<left|{><with|math-font-series|bold|u><rsup|\<asterisk\>(i)><right|}><rsub|i=1><rsup|n-1>> and the associated fluxes <math|<left|{><left|(>A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>(i)><right|}><rsup|n-1><rsub|i=1>>
such that <math|<with|math-font-series|bold|u><rsup|\<asterisk\>(i)>> and
<with|mode|math|<left|(>A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>(i)>>
`lies between' eigenvalues <math|\<lambda\><rsub|i><rsup|\<pm\>>> and
<math|\<lambda\><rsub|i+1><rsup|\<pm\>>> (see figure
<reference|fig:fluxfan>). Now define the unsuperscripted eigenvalues as\
<\equation*>
<tabular|<tformat|<table|<row|<cell|<left|{>\<lambda\><rsub|i><right|}><rsub|i=1><rsup|N>=sort<left|{>\<lambda\><rsub|i><rsup|->:
i\<in\><with|math-font|cal|N><rsup|-><right|}>>>|<row|<cell|>>|<row|<cell|<left|{>\<lambda\><rsub|i><right|}><rsub|i=N+1><rsup|N+Z>=sort<left|{>\<lambda\><rsub|i><rsup|->:
i\<in\><with|math-font|cal|Z><rsup|-><right|}>>>|<row|<cell|>>|<row|<cell|<left|{>\<lambda\><rsub|i><right|}><rsub|i=N+1><rsup|N+Z>=sort<left|{>\<lambda\><rsub|i><rsup|+>:
i\<in\><with|math-font|cal|P><rsup|+><right|}>>>>>>
</equation*>
\;
Define the numerical flux as\
<\equation*>
<left|(>A<with|math-font-series|bold|u><right|)><rsup|#>\<assign\><left|(>A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>(q)>,
</equation*>
where <math|q = argmin<left|{><left|\|>\<lambda\><rsub|i><right|\|><right|}><rsub|i=1><rsup|N>>.\
\;
Assume the following conditions:
<\enumerate>
<item><math|Z\<gtr\>0>
<item>Imposition of the Rankine-Hugoniot condition across eigenvalues:
<\equation*>
\<lambda\><rsub|i><left|(><with|math-font-series|bold|u><rsup|\<asterisk\>(i)>-<with|math-font-series|bold|u><rsup|\<asterisk\>(i-1)><right|)>
= <left|(>A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>(i)>-<left|(>A<with|math-font-series|bold|u><right|)><rsup|\<asterisk\>(i-1)><hspace|1cm>
\<forall\> \ \ i = 1, 2, \<ldots\>n,
</equation*>
where we define <math|<with|math-font-series|bold|u><rsup|\<asterisk\>(0)>=<with|math-font-series|bold|u><rsup|->>
and <math|<with|math-font-series|bold|u><rsup|\<asterisk\>(n)>=<with|math-font-series|bold|u><rsup|+>>.\
<item>Directional biasing:
<\equation*>
(A<with|math-font-series|bold|u>)<rsup|\<asterisk\>(i)>=<left|{><tabular|<tformat|<table|<row|<cell|A<rsup|-><with|math-font-series|bold|u><rsup|\<asterisk\>(i)>>|<cell|>|<cell|if>|<cell|>|<cell|i
= 1, 2, \<ldots\>,N>>|<row|<cell|>|<cell|>|<cell|>|<cell|>|<cell|>>|<row|<cell|A<rsup|+><with|math-font-series|bold|u><rsup|\<asterisk\>(i)>>|<cell|>|<cell|if>|<cell|>|<cell|i
= n-1, n-2, \<ldots\>,n-P>>>>>
</equation*>
<item><math|<with|math-font|cal|W><rsup|+><rsub|<with|math-font|cal|P>>=<with|math-font|cal|W><rsup|-><rsub|<with|math-font|cal|P>>>
and <math|<with|math-font|cal|W><rsub|<with|math-font|cal|N>><rsup|->=<with|math-font|cal|W><rsup|+><rsub|<with|math-font|cal|N>>>
\ (<math|\<Longrightarrow\><with|math-font|cal|W<rsub|Z><rsup|->>=<with|math-font|cal|W><rsub|<with|math-font|cal|Z>><rsup|+>>)
</enumerate>
\;
Then the numerical flux <math|<left|(>A<with|math-font-series|bold|u><right|)><rsup|#>>
exists and is unique. It is given by\
<\equation*>
<tabular|<tformat|<table|<row|<cell|<left|(>A<with|math-font-series|bold|u><right|)><rsup|#>>|<cell|=>|<cell|<big|sum><rsub|i=1><rsup|N>\<nu\><rsub|i>\<lambda\><rsub|i><with|math-font-series|bold|v><rsub|i><rsup|->+<big|sum><rsub|i=n-P+1><rsup|n>\<mu\><rsub|i><rsup|->\<lambda\><rsub|i><with|math-font-series|bold|v><rsub|i><rsup|->>>|<row|<cell|>|<cell|>|<cell|>>|<row|<cell|>|<cell|=>|<cell|<big|sum><rsub|i=1><rsup|N>\<mu\><rsup|+><rsub|i>\<lambda\><rsub|i><with|math-font-series|bold|v><rsub|i><rsup|+>+<big|sum><rsub|i=n-P+1><rsup|n>\<nu\><rsub|i>\<lambda\><rsub|i><with|math-font-series|bold|v><rsub|i><rsup|+>,>>>>>
</equation*>
\;
where the <math|N> unknowns <math|<left|{>\<nu\><rsub|i><right|}><rsub|i=1><rsup|N>>
satisfy the rank <math|N> system\
<\equation*> <left|(><tabular|<tformat|<table|<row|<cell|\<lambda\><rsub|1><with|math-font-series|bold|v><rsub|1><rsup|->>|<cell|\<lambda\><rsub|2><with|math-font-series|bold|v><rsub|2><rsup|->>|<cell|>|<cell|\<cdots\>>|<cell|>|<cell|\<lambda\><rsub|N><with|math-font-series|bold|v><rsub|N><rsup|->>>>>><right|)><left|(><tabular|<tformat|<table|<row|<cell|\<nu\><rsub|1>>>|<row|<cell|\<nu\><rsub|2>>>|<row|<cell|\<vdots\>>>|<row|<cell|\<nu\><rsub|N>>>>>><right|)>=<big|sum><rsub|i=1><rsup|N>\<mu\><rsup|+><rsub|i>\<lambda\><rsub|i><with|math-font-series|bold|v><rsub|i><rsup|+>,
</equation*>
and the <math|P> unknowns <math|<left|{>\<nu\><rsub|i><right|}><rsub|i=n-P+1><rsup|n>>
satisfy the rank <math|P> system\
<\equation*>
<left|(><tabular|<tformat|<table|<row|<cell|\<lambda\><rsub|n-P+1><with|math-font-series|bold|v><rsub|n-P+1><rsup|+>>|<cell|\<lambda\><rsub|n-P+2><with|math-font-series|bold|v><rsub|n-P+2><rsup|+>>|<cell|>|<cell|\<cdots\>>|<cell|>|<cell|\<lambda\><rsub|n><with|math-font-series|bold|v><rsub|n><rsup|+>>>>>><right|)><left|(><tabular|<tformat|<table|<row|<cell|\<nu\><rsub|n-P+1>>>|<row|<cell|\<nu\><rsub|n-P+2>>>|<row|<cell|\<vdots\>>>|<row|<cell|\<nu\><rsub|n>>>>>><right|)>=<big|sum><rsub|i=n-P+1><rsup|n>\<mu\><rsup|-><rsub|i>\<lambda\><rsub|i><with|math-font-series|bold|v><rsub|i><rsup|->.
</equation*>
</theorem>
Note that this theorem even takes into account degeneracies in any
eigenvalue.
<subsection|Deriving an Upwind Flux for the Wave Equation>
We're dealing with
<\equation*>
u<rsub|t t>=c<rsup|2>(x)\<Delta\>u,
</equation*>
which we rewrite into conservation form:
<\eqnarray*>
<tformat|<table|<row|<cell|u<rsub|t>-c(x)\<nabla\>\<cdot\>\<b-v\>>|<cell|=>|<cell|0>>|<row|<cell|\<b-v\><rsub|t>-c(x)\<nabla\>u>|<cell|=>|<cell|0>>>>
</eqnarray*>
and from there into matrix form:
<\eqnarray*>
<tformat|<table|<row|<cell|\<b-w\><rsub|t>+c(x)<left|[><wide*|<matrix|<tformat|<table|<row|<cell|0>|<cell|-1>|<cell|0>>|<row|<cell|-1>|<cell|>|<cell|>>|<row|<cell|0>|<cell|>|<cell|>>>>>|\<wide-underbrace\>><rsub|A<rsub|x>\<assign\>>\<b-w\><rsub|x>+<wide*|<matrix|<tformat|<table|<row|<cell|0>|<cell|0>|<cell|-1>>|<row|<cell|0>|<cell|>|<cell|>>|<row|<cell|-1>|<cell|>|<cell|>>>>>|\<wide-underbrace\>><rsub|A<rsub|y>\<assign\>>\<b-w\><rsub|y><right|]>>|<cell|=>|<cell|\<b-0\>>>>>
</eqnarray*>
where we set <math|\<b-w\>\<assign\>(u,\<b-v\>)>.\
(see <with|font-family|tt|doc/maxima/wave.mac> and
<with|font-family|tt|doc/maxima/myhelpers.mac> for Maxima scripts for code
to find the upwind flux.)
We find:
<\eqnarray*>
<tformat|<table|<row|<cell|(\<Pi\>*u)<rsup|\<ast\>>>|<cell|=>|<cell|\<b-n\>\<cdot\><average|c\<b-v\>>+<frac|[c*u]<rsup|->-[c*u]<rsup|+>|2>>>|<row|<cell|(\<Pi\>\<b-v\>)<rsup|\<ast\>>>|<cell|=>|<cell|\<b-n\><average|c*u>+<frac|1|2>(\<b-n\>\<otimes\>\<b-n\>)\<cdot\>[c<rsup|->\<b-v\><rsup|->-c<rsup|+>\<b-v\><rsup|+>].>>>>
</eqnarray*>
<subsection|Deriving an Upwind Flux for the 2D Advection Equation>
We're dealing with
<\equation*>
u<rsub|t>=\<b-v\>\<cdot\>\<nabla\>u,
</equation*>
which we rewrite into matrix conservation form as
<\equation*>
u<rsub|t>=v<rsub|1>u<rsub|x>+v<rsub|2>u<rsub|y>.
</equation*>
We find <math|A(\<b-n\>)=\<b-n\>\<cdot\>\<b-v\>>, and therefore
<\equation*>
(A(\<b-n\>)u)<rsup|\<ast\>>=<choice|<tformat|<table|<row|<cell|\<b-n\>\<cdot\>\<b-v\>u<rsup|->>|<cell|<with|mode|text|if
<math|\<b-n\>\<cdot\>\<b-v\>\<geqslant\>0>
(outflow)>,>>|<row|<cell|\<b-n\>\<cdot\>\<b-v\>u<rsup|+>>|<cell|<with|mode|text|if
<math|\<b-n\>\<cdot\>\<b-v\>\<less\>0> (inflow)>.>>>>>
</equation*>
<section|Tim's funky curvilinear element stuff>
We're discretizing <with|mode|math|u<rsub|t>+v\<cdot\>\<nabla\>u=0>. Recall
<\eqnarray*>
<tformat|<table|<row|<cell|\<nabla\>\<cdot\>(v*u*\<varphi\>)>|<cell|=>|<cell|\<nabla\>\<cdot\>(v*u)\<varphi\>+v*u\<cdot\>\<nabla\>\<varphi\>=(v\<cdot\>\<nabla\>*u)\<varphi\>+v*u\<cdot\>\<nabla\>\<varphi\>.>>>>
</eqnarray*>
Using the summation convention, we find
<\eqnarray*>
<tformat|<table|<row|<cell|0>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|T<rsub|k>>(v\<cdot\>\<nabla\>u)\<varphi\>>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>v*u\<cdot\>\<nabla\>\<varphi\>+<big|int><rsub|T<rsub|k>>\<nabla\>\<cdot\>(v*u*\<varphi\>)>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>v*u\<cdot\>\<nabla\>\<varphi\>+<big|int><rsub|\<partial\>T<rsub|k>>v*u\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|\<approx\>>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>v*u\<cdot\>\<nabla\>\<varphi\>+<big|int><rsub|\<partial\>T<rsub|k>>(v*u)<rsup|\<ast\>>\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|\<approx\>>|<cell|<big|int><rsub|T<rsub|k>>\<partial\><rsub|t>u<rsub|i>l<rsub|i>l<rsub|j>-<big|int><rsub|T<rsub|k>><matrix|<tformat|<table|<row|<cell|a<rsub|1>u<rsub|i>l<rsub|i>>>|<row|<cell|a<rsub|2>u<rsub|i>l<rsub|i>>>>>>\<cdot\><matrix|<tformat|<table|<row|<cell|\<partial\><rsub|x<rsub|2>>l<rsub|j>>>|<row|<cell|\<partial\><rsub|x<rsub|2>>l<rsub|j>>>>>>+<big|sum><rsub|F\<subset\>\<partial\>T<rsub|k>><big|int><rsub|F><matrix|<tformat|<table|<row|<cell|(v<rsub|1>u<rsub|i>)<rsup|\<ast\>>l<rsub|i><rsup|F>l<rsub|j><rsup|F>>>|<row|<cell|(v<rsub|2>u<rsub|i>)<rsup|\<ast\>>l<rsub|i><rsup|F>l<rsub|j><rsup|F>>>>>>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|>>|<row|<cell|>|<cell|=>|<cell|>>>>
</eqnarray*>
<\eqnarray*>
<tformat|<table|<row|<cell|0>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>-<big|int><rsub|D<rsub|k>>a*u\<cdot\>\<nabla\>\<varphi\>+<big|int><rsub|\<partial\>T<rsub|k>>(v*u)<rsup|\<ast\>>\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|D<rsub|k>>(v\<cdot\>\<nabla\>u)\<varphi\>-<big|int><rsub|\<partial\>T<rsub|k>>(v*u-{v*u})\<varphi\>\<cdot\>n>>|<row|<cell|>|<cell|=>|<cell|<big|int><rsub|T<rsub|k>>u<rsub|t>\<varphi\>+<big|int><rsub|D<rsub|k>>(v\<cdot\>\<nabla\>u)\<varphi\>-<big|int><rsub|\<partial\>T<rsub|k>>v<matrix|<tformat|<table|<row|<cell|<frac|1|2>>>|<row|<cell|-<frac|1|2>>>>>>\<cdot\><matrix|<tformat|<table|<row|<cell|u<rsup|->>>|<row|<cell|u<rsup|+>>>>>>\<varphi\>\<cdot\>n>>>>
</eqnarray*>
</body>
<\initial>
<\collection>
<associate|info-flag|detailed>
<associate|page-orientation|portrait>
<associate|page-type|letter>
<associate|par-first|0>
</collection>
</initial>
<\references>
<\collection>
<associate|auto-1|<tuple|1|1>>
<associate|auto-10|<tuple|5.3|7>>
<associate|auto-11|<tuple|5.3.1|7>>
<associate|auto-12|<tuple|5.3.2|8>>
<associate|auto-13|<tuple|5.4|9>>
<associate|auto-14|<tuple|6|9>>
<associate|auto-15|<tuple|7|10>>
<associate|auto-16|<tuple|7.1|10>>
<associate|auto-17|<tuple|1|11>>
<associate|auto-18|<tuple|7.2|11>>
<associate|auto-19|<tuple|7.2.1|14>>
<associate|auto-2|<tuple|2|1>>
<associate|auto-20|<tuple|7.2.2|14>>
<associate|auto-21|<tuple|7.3|18>>
<associate|auto-22|<tuple|7.3.1|?>>
<associate|auto-23|<tuple|7.3.2|?>>
<associate|auto-24|<tuple|7.4|?>>
<associate|auto-25|<tuple|7.5|?>>
<associate|auto-26|<tuple|8|?>>
<associate|auto-27|<tuple|8.0.1|?>>
<associate|auto-28|<tuple|8|?>>
<associate|auto-29|<tuple|8|?>>
<associate|auto-3|<tuple|3|3>>
<associate|auto-4|<tuple|4|4>>
<associate|auto-5|<tuple|5|6>>
<associate|auto-6|<tuple|5.1|6>>
<associate|auto-7|<tuple|5.2|7>>
<associate|auto-8|<tuple|5.2.1|7>>
<associate|auto-9|<tuple|5.2.2|7>>
<associate|auto.1-1|<tuple|1|?|#1>>
<associate|auto.2-1|<tuple|2|?|#2>>
<associate|auto.3-1|<tuple|3|?|#3>>
<associate|auto.4-1|<tuple|4|?|#4>>
<associate|auto.5-1|<tuple|5|?|#5>>
<associate|auto.5-2|<tuple|5.1|?|#5>>
<associate|auto.5-3|<tuple|5.2|?|#5>>
<associate|auto.5-4|<tuple|5.2.1|?|#5>>
<associate|auto.5-5|<tuple|5.2.2|?|#5>>
<associate|auto.5-6|<tuple|5.3|?|#5>>
<associate|auto.5-7|<tuple|5.3.1|?|#5>>
<associate|auto.5-8|<tuple|5.3.2|?|#5>> <associate|auto.5-9|<tuple|5.4|?|#5>>
<associate|auto.6-1|<tuple|6|?|#6>>
<associate|auto.7-1|<tuple|7|?|#7>>
<associate|auto.7-10|<tuple|7.4|?|#7>>
<associate|auto.7-11|<tuple|7.5|?|#7>>
<associate|auto.7-12|<tuple|7.6|?|#7>>
<associate|auto.7-13|<tuple|7.7|?|#7>>
<associate|auto.7-2|<tuple|7.1|?|#7>>
<associate|auto.7-3|<tuple|1|?|#7>>
<associate|auto.7-4|<tuple|7.2|?|#7>>
<associate|auto.7-5|<tuple|7.2.1|?|#7>>
<associate|auto.7-6|<tuple|7.2.2|?|#7>>
<associate|auto.7-7|<tuple|7.3|?|#7>>
<associate|auto.7-8|<tuple|7.3.1|?|#7>>
<associate|auto.7-9|<tuple|7.3.2|?|#7>>
<associate|auto.8-1|<tuple|8|?|#8>>
<associate|auto.9-1|<tuple|9|?|#9>>
<associate|eq:rh-all|<tuple|1|?>>
<associate|eq:rh-aux|<tuple|2|?>>
<associate|fig:fluxfan|<tuple|1|11>>
</collection>
</references>
<\auxiliary>
<\collection>
<\associate|figure>
<tuple|normal|<label|fig:fluxfan>Space-Time Diagram of a flux
fan.|<pageref|auto-17>>
</associate>
<\associate|toc>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|1<space|2spc>Mapping
from 2D equilateral to Unit Coordinates>
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-1><vspace|0.5fn>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|2<space|2spc>Mapping
from 3D equilateral to Unit Coordinates>
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-2><vspace|0.5fn>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|3<space|2spc>Derivatives
of the 2D Basis functions> <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-3><vspace|0.5fn>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|4<space|2spc>Derivatives
of the 3D Basis functions> <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-4><vspace|0.5fn>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|5<space|2spc>DG
Schemes> <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-5><vspace|0.5fn>
<with|par-left|<quote|1.5fn>|5.1<space|2spc>Notation
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-6>>
<with|par-left|<quote|1.5fn>|5.2<space|2spc>Advection Equation
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-7>>
<with|par-left|<quote|3fn>|5.2.1<space|2spc>Weak DG
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-8>>
<with|par-left|<quote|3fn>|5.2.2<space|2spc>Strong DG
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-9>>
<with|par-left|<quote|1.5fn>|5.3<space|2spc>Wave Equation
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>> <no-break><pageref|auto-10>>
<with|par-left|<quote|3fn>|5.3.1<space|2spc>Weak DG
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-11>>
<with|par-left|<quote|3fn>|5.3.2<space|2spc>Strong DG
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-12>>
<with|par-left|<quote|1.5fn>|5.4<space|2spc>Heat Equation
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-13>>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|6<space|2spc>Quadrature
Rules> <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-14><vspace|0.5fn>
<vspace*|1fn><with|font-series|<quote|bold>|math-font-series|<quote|bold>|7<space|2spc>Upwind
Fluxes> <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-15><vspace|0.5fn>
<with|par-left|<quote|1.5fn>|7.1<space|2spc>DIY Upwind Fluxes/Riemann
Solvers <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-16>>
<with|par-left|<quote|1.5fn>|7.2<space|2spc>Comments by Akil
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-18>>
<with|par-left|<quote|3fn>|7.2.1<space|2spc>Zero eigenvalues
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-19>>
<with|par-left|<quote|3fn>|7.2.2<space|2spc>Degenerate eigenvalues
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-20>>
<with|par-left|<quote|1.5fn>|7.3<space|2spc>More random comments by
Akil <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-21>>
<with|par-left|<quote|3fn>|7.3.1<space|2spc>Stoopid upwinding
<datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-22>>
<with|par-left|<quote|3fn>|7.3.2<space|2spc>R-H Conditions with
Directional Bias <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-23>>
<with|par-left|<quote|1.5fn>|7.4<space|2spc>Deriving an Upwind Flux for
the Wave Equation <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-24>>
<with|par-left|<quote|1.5fn>|7.5<space|2spc>Deriving an Upwind Flux for
the 2D Advection Equation <datoms|<macro|x|<repeat|<arg|x>|<with|font-series|medium|<with|font-size|1|<space|0.2fn>.<space|0.2fn>>>>>|<htab|5mm>>
<no-break><pageref|auto-25>>
</associate>
</collection>
</auxiliary>