Excercise sheet 10 -- solutions

Matrices

1. Define matrices
	A=(({{1, 2}, {3, 4}})), B=(({{0, 1}, {1, 1}})),
and the vectors 
	X=(({{1}, {2}})), Y=(({{3}, {4}})).
Do the following
(a) Calculate the expressions A B,  B A ,  A X,  
      X·Y and A B - B A
(b) Find the determinant of A, det A.
(c) Solve the equation A Z = Xfor Z. Check the result by substituting Z back into the equation
(d) Find the inverse (A^(-1)) and calculate (A^(-1)) A.
(e) Find the transposed matrix (A^T).
(f) Find the eigenvalues and eigenvectors of A.
(g) Using Mathematicas help system, can you find a command that calculates the cross-product? Use it to calculate X×Y.

2. Vector calculus Mathematica commands can be found in the package Calculus`VectorAnalysis`. Can you find a way to calculate the gradient ∇f(x, y, z), divergence \(TraditionalForm\`∇\(\(·\)\(\(G\&_\)(x, y, z)\)\)\) and curl \(TraditionalForm\`∇\(\(×\)\(\(G\&_\)(x, y, z)\)\)\) for the scalar function f(x, y, z) = 1/(Sqrt[(x^2) + (y^2) + (z^2)]) and vector function \(TraditionalForm\`\(G\&_\)(x, y, z) = \((x\^2 + y, y + z, z\ x)\)\)?

Solutions

1

A = {{1, 2}, {3, 4}}; B = {{0, 1}, {1, 1}}; X = {1, 2}; 
    Y = {3, 4};

Print["AB = ", 
      A . B // MatrixForm]

    Print["BA = ", B . A // MatrixForm]

    Print["AX = ", A . X // MatrixForm]

    Print["X·Y = ", X . Y // MatrixForm]

    Print["AB-BA = ", 
      A . B - B . A // 
        MatrixForm, ". This is the commutator of A and B."]

Det[A]

Z = 
        LinearSolve[A, X];

    Print["Z = ", Z // MatrixForm]

    Print["AZ = ", A . Z // MatrixForm, ", X = ", 
      X // MatrixForm]

\(\(\( (*d*) \)\(
\)\(Print[\*""(A^(-1)) = "", Inverse[A] // MatrixForm]

    Print[\*""(A^(-1))A = "", 
      Inverse[A] . A // MatrixForm, \*"", (AA^(-1)) = "", 
      A . Inverse[A] // MatrixForm]\)\)\)

\(\(\( (*e*) \)\(
\)\(Print[\*""\!\(A\^T\) = "", 
      Transpose[A] // MatrixForm]\)\)\)

Print["The eigenvalues of A are \
", Eigenvalues[A][[1]], " and ", Eigenvalues[
          A][[2]], "."]

    Print["The corresponding eigenvectors are ", Eigenvectors[
            A][[1]] // 
        MatrixForm, " and ", Eigenvectors[A][[2]] // 
        MatrixForm, "."]

In (g), we'll assume that the third component of X and Y is zero (in order to get the familiar cross product).

X = Append[X, 0]; 
    Y = Append[Y, 0];

    Print["X\[Cross]Y = ", Cross[X, Y] // MatrixForm]

2

<< Calculus`VectorAnalysis`

The syntax of the vectors analysis commands can be found in the help browser. First, we can set the default coordinate system:

SetCoordinates[Cartesian[x, y, z]];

f[x_, y_, z_] = 
        (1/(Sqrt[(x^2) + (y^2) + (z^2)]));
G[x_, y_, z_] = {(x^2) + y, y + z, z x};

The nabla symbol is typed Esc-del-Esc.

Print["∇f(x,y,z) = ", 
      Grad[f[x, y, z]]]
Print["∇·G(x,y,z) = ", 
      Div[G[x, y, z]]]
Print["∇\[Cross]G(x,y,z) = ", Curl[G[x, y, z]]]

Schrödinger equation in matrix form

In the lecture notes, the Schrödinger equation was solved for the harmonic potential. Modify the potential to include the term ϵ (r^4), that is, use V(r) = (r^2) + ϵ (r^4)  Set ϵ = 0.1. How do the energies change compared to the original, harmonic potential? How about the eigenfunctions?

Solution

The equation in question is y'' (r) + V(r) y(r) = E y(r). The idea behind the matrix way of solving a differential equation is to consider the left hand side of this equation as an operator acting on the function y, i.e. L[y] = E y, and then transform the operator into matrix form with the method of finite differences. Then we have a matrix eigenvalue equation, denoted by M Y = E Y, which is easy to solve using a computer mathematics program. The eigenvectors of M, i.e. the vectors Y, are vectors with the values of the function y as their elements, i.e. ((Y = (y((Subscript[r, 0])), y((Subscript[r, 1])),  ... , y((Subscript[r, n]))))^T). The second derivative of y is transformed as
	\(TraditionalForm\`y'' \((r\_i)\) = \(\(y(r\_\(i + 1\)) - 2 \( y(r\_i)\) + y(r\_\(i - 1\))\)\/h\^2 = \(y\_\(i + 1\) - 2  y\_i + y\_\(i - 1\)\)\/h\^2\)\),
so Schrödinger's equation in finite difference form looks like
	-(((Subscript[y, i + 1]) - 2  (Subscript[y, i]) + 
                  (Subscript[y, i - 1]))/(h^2)) + (Subscript[V, i]) (Subscript[y, i]) = E (Subscript[y, i])
-(1/(h^2)) (Subscript[y, i - 1]) + ((2/(h^2)) + (Subscript[V, i])) 
            (Subscript[y, i]) - (1/(h^2)) (Subscript[y, i + 1]) = E (Subscript[y, i])
The idea is to form this equation at each of the points {(Subscript[r, i])}. Then we have a set of n linear equations for the variables (Subscript[y, i]), which can be written as M Y  =  E Y, and the elements of the matrix M are \(TraditionalForm\`\(m\_\(i, j\) = \(\*"""the\ coefficient\ of\ y\_j\ in\ the\ i : th\ equation\)\*"""\)\). For example (Subscript[m, i, i + 1]) = -(1/(h^2)), (Subscript[m, i, i]) = (2/(h^2)) + (Subscript[V, i]). Most of the elements of M are zero. In fact, only the diagonal and the adjacent bands contain non-zero elements:
M=(({{(Subscript[m, 11]), (Subscript[m, 12]), 0, 0, 0}, {(Subscript[m, 21]), (Subscript[m, 22]), (Subscript[m, 23]), 0, 0}, {0, (Subscript[m, 32]), (Subscript[m, 33]), (Subscript[m, 34]), 0}, {0, 0, (Subscript[m, 43]), (Subscript[m, 44]), (Subscript[m, 45])}, {0, 0, 0, (Subscript[m, 54]), (Subscript[m, 55])}})).

First, define the interval [(Subscript[r, a]), (Subscript[r, b])] in which to search for the solution, and the spacing of the points {(Subscript[r, 0]) = (Subscript[r, a]), (Subscript[r, 1]) = (Subscript[r, a]) + h, 
        (Subscript[r, 2]) = (Subscript[r, a]) + 2  h,  ... , (Subscript[r, n]) = (Subscript[r, b])} on this interval:

h = 0.05;
ra = 0.0;
rb = 20.0;
n = (((Subscript[r, b]) - (Subscript[r, a]))/h) + 1;

The potential is of the form V(r) = (r^2) + ϵ (r^4), with ϵ = 0.1. Calculate its values in points {(Subscript[r, i])}:

ϵ = 0.1;
V[r_] := (r^2) + ϵ (r^4);
vtab = Table[V[r], {r, ra, rb, h}];

Then, construct the matrix M. The only non-zero elements for the ith row are (Subscript[m, i, i - 1]), (Subscript[m, i, i]) and (Subscript[m, i, i + 1]).

M = 
        Table[Switch[i - j, -1, -(1/(h^2)), 1, -(1/(h^2)), 0, 
            (2/(h^2)) + vtab[[i]], _, 0], {i, n}, {j, n}];

Now, the solutions of the equation M Y = E Y, the eigenvalues E and the corresponding eigenvectors Y can be calculated:

eigvec = Eigenvectors[M];
eigval = Eigenvalues[M];

Then, we can plot the eigenfunctions, which can be obtained by means of interpolation.

np = 5;
Plot[
        Evaluate[
          Table[Interpolation[eigvec[[n - i]]][(x - ra)/h + 
                1], {i, 0, np}]], {x, ra, 5}, PlotRange -> All, 
        PlotStyle -> Table[Hue[i/np], {i, 0, np}], 
        ImageSize -> 700];