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!