10.5 There and back again - an exploration of the liouville  (Page 9/10)

In computing the results which are reported starting in "Numerical Results" , various techniques were used. Some techniques were successful and other techniques failed to give the desired or expected results, which led to alternate techniques being tried. A short description of our methods, including a description of spectral discretization, is given in the next section.

Numerical methods

When we began gathering the data included here in our results, we started by programming the $q\left(t\right)$ in [link] from [link] into MATLAB. Our next goal was to find the $\rho \left(x\right)$ from these potential functions. We started by assuming $x\in [0,1]$ and used a uniform grid to solve the ODE from the Liouville transformation [link] with forward Euler and the change of variables for $t$ [link] with composite trapezoid quadrature. Briefly, forward Euler gives a numerical approximation to the solution of the initial value problem

through the formula

where $h=b_{n+1}-b_n$ [link] . In order to facilitate finding the eigenvalues, we attempted to incorporate a Chebyshev grid for $x$ into the uniform grid. This is when we became aware something was not working properly. In an attempt to find the problem, we output the value of $t$ at each step, and we realized for some choices of ${\mu }_1$ and initial conditions, when $x$ reached 1, $t$ might be much less than 1 or even greater than 1. We had two choices at this point, we could choose values of ${\mu }_1$ and initial conditions and find the necessary value of $x$ which would result in $t$ being within some tolerance of 1 or we could solve the ODE in $t$ instead of $x$ . This second method is how we chose to continue our exploration.

With a uniform grid over $t\in [0,1]$ and given ${\mu }_1$ , $\rho \left(0\right)$ , and ${\rho }^{\text{'}}\left(0\right)$ we used forward Euler to solve the ODE [link] and then implemented the change of variables. This resulted in strings of various lengths depending on the choice of ${\mu }_1$ , $\rho \left(0\right)$ , and ${\rho }^{\text{'}}\left(0\right)$ . These are illustrated in [link] .

We made the choice to restrict to strings of unit length and $\rho \left(0\right)=1$ . To facilitate this choice, we included a bisection solver to find the necessary value of ${\rho }^{\text{'}}\left(0\right)$ for a given ${\mu }_1$ that results in a string of length 1. After this determination is made, the resulting graphs are illustrated in [link] .

For a given ${\mu }_1$ by changing the choice of $\rho \left(0\right)$ , with the necessary value of ${\rho }^{\text{'}}\left(0\right)$ being determined to result in a unit string, the isospectral graphs, graphs with the same Dirichlet spectrum, are produced (seen in [link] ).

To compute the eigenvalues and eigenfunctions, we wanted to use spectral discretization, but because our method solves the ODE in $t$ , we could not choose at which $x$ values we would have values for $\rho \left(x\right)$ . To overcome this, we used the MATLAB function "interp1" and used a cubic spline interpolation to find $\rho$ at the Chebyshev grid x values from the values we did have from solving the ODE in $t$ . Then, after adjusting the differential operator matrix, $D$ , and $\rho$ matrix, $B$ , for the Dirichlet boundary conditions, we used the MATLAB function "eig" to solve for $M$ , the diagonal matrix of generalized eigenvalues, and $V$ , a full matrix of corresponding eigenvectors, in the generalized eigenvalue problem, $D^2*V=M*B*V$ (from $-y^{\text{'}\text{'}}\left(x\right)=\lambda \rho \left(x\right)y\left(x\right)$ ). The eigenvalues were sorted, and the eigenvectors were sorted correspondingly, to produce the figures in this paper.