Solving system of ordinary differential equations where coefficients are matrices

Question

I am attempting to use NDSolve to solve a system of ordinary differential equations defined by matrices. Here's the setup of my problem and what I have tried so far.

I'm interested in applying the Galerkin Method to solve Burger's 1-D viscous equation using Chebyshev polynomials as test functions. My test functions are therefore defined as:

ψ[i_, x_] := ChebyshevT[i, x]

And I am looking for a solution that is of the form:

$$ u_{a}(x,t) = \sum_{j=0}^{N}c_{j}(t)T_{j}(x), $$

where $T_{j}(x)$ are the Chebyshev polynomials I'm using as the test function. After finding the residual and projecting it onto the basis functions $(R,T_{k})=0$, I get my first equation to solve:

$$ \underline{\underline{M}}\underline{\dot{c}}+\underline{\underline{B}}\underline{c}+\underline{\underline{C}}\underline{c}=0. $$

My second equation is found by solving for $\underline{\dot{c}}$:

$$ \underline{\dot{c}}=-\underline{\underline{M}}^{-1}\underline{\underline{B}}\underline{c}-\underline{\underline{M}}^{-1}\underline{\underline{C}}\underline{c}. $$

The matrices $\underline{\underline{M}}, \underline{\underline{B}}, \: \& \: \underline{\underline{C}}$ are the inner products of the residual projected onto my test function, and I have programmed them as:

$\underline{\underline{M}}$:

mjk[j_, k_] := mjk[j, k] = Integrate[(ψ[j, x] ψ[k, x])/Sqrt[1 - x^2], {x, -1, 1}]
m2jk[n_] := m2jk[n] = Table[Table[mjk[j, k], {j, 1, n}], {k, 1, n}]

$\underline{\underline{B}}$:

bjk[j_, k_, i_] := bjk[j, k, i] = Integrate[ψ[j, x] D[ψ[k, x], {x, 1}] ψ[i, x]/Sqrt[1 - x^2], {x, -1, 1}]
b2jk[n_] := b2jk[n] = Table[Table[bjk[j, k, n], {j, 1, n}], {k, 1, n}]

$\underline{\underline{C}}$:

cjk[j_, k_] := cjk[j, k] = Integrate[(-1/RE) D[ψ[j, x], {x, 1}] D[ψ[k, x], {x, 1}]/Sqrt[1 - x^2], {x, -1, 1}]
c2jk[n_] := c2jk[n] = Table[Table[cjk[j, k], {j, 1, n}], {k, 1, n}]

Finally, I try to solve the system of ODEs for the coefficients $c$ in the following way:

odes1[n_] := Table[
   Evaluate[m2jk[n].Table[D[γ[k, t], t], {k, 1, n}]][[jj]]
   + Evaluate[b2jk[n].Table[γ[k, t], {k, 1, n}]][[jj]]
   + Evaluate[c2jk[n].Table[γ[k, t], {k, 1, n}]][[jj]] == 0,
 {jj, 1, n}]

odes2[n_] :=Table[
   Evaluate[Table[D[γ[k, t], t], {k, 1, n}]][[jj]]
     == Evaluate[ccg[n] bb[n].Table[γ[k, t], {k, 1, n}]][[jj]]
        + Evaluate[cc[n].Table[γ[k, t], {k, 1, n}]][[jj]],
 {jj, 1, n}]

c[n_] := Table[γ[k, 0] == 0, {k, 1, n}]

solns[n_] := Part[NDSolve[{odes2[n], c[n]}, Table[γ[jj, t], {jj, 1, n}], {t, 0, 1}], 1]

solutions[n_] := Table[γ[kk, t], {kk, 0, n}] /. solns[n]

Plot[Evaluate[solutions[3]], {t, 0, 1}, PlotRange -> All]

But as can be seen by the above plot, the result that I am getting is zero (which I know to be incorrect). I don't think that I made an error in my definition of the matrices $\underline{\underline{M}}, \underline{\underline{B}}, \: \& \: \underline{\underline{C}}$, but I am unsure of my utilization of NDSolve. I hope that I have given enough background to the problem that I am attempting to solve (for anyone curious, this is example 1.2.5 in "Computational Galerkin Methods" by C.A.J. Fletcher). Any tips as to where I am going wrong in my attempt to use NDSolve to solve the system of ODEs for the unknown $c$'s would be greatly appreciated. Thank you in advance for your time!

Answer 1

Thanks to @Ulrich Neumann for pointing out a missing parameter definition; it caused me to look at this again and see quite a few problems with my solution methodology! Therefore I hope that answering my own question will be allowed in the hopes that it gives some insight to someone else.

To start, the system of equations formed by substituting the approximate solution into the equation governing Burgers' 1D problem and projecting the resultant residual onto the Chebyshev polynomials (i.e. the basis functions) stated in the original problem statement was incomplete. Following along with section 1.2.5. in Fletcher's "Computational Galerkin Methods" text, I tried to replicate the code "BURG1".

To solve the system of equations, the following matrices were constructed:

mjk[j_, k_] := 
mjk[j, k] = Integrate[\[Psi][j, x] \[Psi][k, x], {x, -1, 1}]
m2jk[n_] := m2jk[n] = Table[Table[mjk[j, k], {j, 1, n}], {k, 1, n}]

oj[j_] := oj[j] = Integrate[(-(1 - x)/4) \[Psi][j, x], {x, -1, 1}]
o2j[n_] := o2j[n] = Table[oj[j], {j, 1, n}]

njk[j_, k_] := njk[j, k] = 
Integrate[((1 - x)/2) D[\[Psi][j, x], {x, 1}] \[Psi][k, x], {x, -1, 1}]
n2jk[n_] := n2jk[n] = Table[Table[njk[j, k], {j, 1, n}], {k, 1, n}]

qjk[j_, k_, Re_] := qjk[j, k, Re] = 
Integrate[(-1/Re) \[Psi][j, x] \[Psi][k, x], {x, -1, 1}]
q2jk[n_, Re_] := 
q2jk[n, Re] = Table[Table[qjk[j, k, Re], {j, 1, n}], {k, 1, n}]

bjk[j_, k_] := bjk[j, k] = 
Integrate[\[Psi][j, x] D[\[Psi][j, x], {x, 1}] \[Psi][k, x], {x, -1, 1}]
b2jk[n_] := b2jk[n] = Table[Table[bjk[j, k], {j, 1, n}], {k, 1, n}]

cjk[j_, k_, Re_] := cjk[j, k, Re] = 
Integrate[(-1/Re) D[\[Psi][j, x], {x, 2}] \[Psi][k, x], {x, -1, 1}]
c2jk[n_, Re_] := 
c2jk[n, Re] = Table[Table[cjk[j, k, Re], {j, 1, n}], {k, 1, n}]

The next bit was to make the ODE equation a bit easier to follow...

mmi[n_] := Inverse[m2jk[n]]
oo[n_] := -mmi[n].o2j[n]
nn[n_] := -mmi[n].n2jk[n]
qq[n_, Re_] := -mmi[n].q2jk[n, Re]
bb[n_] := -mmi[n].b2jk[n]
cc[n_, Re_] := -mmi[n].c2jk[n, Re]

Finally, the second equation solving for $ c_{j}(t) $ can be formed as:

odes[n_, Re_] :=Table[
Evaluate[Table[D[\[Gamma][k, t], t], {k, 1, n}]][[jj]]
== Evaluate[oo[n]][[jj]]
+ Evaluate[nn[n].Table[\[Gamma][k, t], {k, 1, n}]][[jj]]
+ Evaluate[qq[n, Re].Table[\[Gamma][k, t], {k, 1, n}]][[jj]]
+ Evaluate[bb[n].Table[\[Gamma][k, t], {k, 1, n}]][[jj]]
+ Evaluate[cc[n, Re].Table[\[Gamma][k, t], {k, 1, n}]][[jj]], {jj, 1, n}]

Solving for the unknown coefficients $c_{j}(t)$ using an initial solution as:

cjt[n_] := Table[\[Gamma][i, 0] == 
Integrate[(u0[x] - ua0[x]) \[Psi][i, x], {x, -1, 1}], {i, 1, n}]

Solving for the approximate solution can be accomplished by:

solns[n_, Re_] :=Part[NDSolve[{odes[n, Re], cjt[n]},
Table[\[Gamma][jj, t], {jj, 1, n}], {t, 0, 1}], 1]

solutions[n_, Re_] := Table[\[Gamma][kk, t], {kk, 1, n}] /. solns[n, Re]

Subscript[u, a][x_, n_, Re_] := ua0[x] +
Sum[Evaluate[Table[\[Gamma][kk, t], {kk, 1, n}] /. solns[n, Re]]
[[jj]] \[Psi][jj, x], {jj, 1, n}]

u[x_, tt_, n_, Re_] := Subscript[u, a][x, n, Re] /. t -> tt

Comparing with the analytical solution to Burgers' 1D solution at a Reynolds number of 10 for N = 4 terms, the approximate solution follows the exact solution, but there are still significant differences...

It is quite easy to track the propagation of the "shock" with time as shown below...

For higher terms (N = 7, 9, and 11), comparison with the results published in Table 1.8 of Fletcher's "Computational Galerkin Methods" shows quite a bit of discrepency.

Additionally, the "shock" still doesn't appear quite right to me; it should be a bit more sudden and not so dissipative. So this is maybe not a full and proper answer, but I think it helps to answer the original problem that I asked, which is my reason for posting.