Excercise 10

1. Matrices

1.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.

1.2 Advanced - try this only if you are interested and have time. 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.1

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

A . B

    B . A

    A . X

    X . Y

    A . B - B . A
  

    Det[A]
  

    LinearSolve[A, X]

    A . % == X
  

    Inverse[A]

    % . A
  

    Transpose[A]
  

    Eigenvalues[A]

    Eigenvectors[A]
  

    Cross[{1, 0, 0}, {0, 1, 0}]

1.2

Clear[f, G, x, y, z]
<< Calculus`VectorAnalysis`
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}

Grad[f[x, y, z]]
Div[G[x, y, z]]
Curl[G[x, y, z]]

2. Schrödinger equation in matrix form

2.1 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

First of all, here are the commands to find the three lowest energies for the harmonic potential. This piece of code is taken from the lecture notes, and slightly improved.

h   =  0.05;  
alku  =  0; 
loppu =  8; 
nn := (loppu - alku)/h + 1;
Clear[x, r, v, i, j, d]
v = Table[r^2, {r, alku, loppu, h}];
d = -1/h^2
mat = 
        Table[Switch[i - j, -1, d, 0, v[[i]] -  2 d, 1, d, _, 
            0], {i, nn}, {j, nn}];
  
eigen = Reverse[Eigenvalues[mat, -10]]   
eigvec = 
        Reverse[Eigenvectors[
            mat, -10]];

ListPlot[eigen]
lp1 = 
      ListPlot[eigvec[[1]], PlotJoined -> True, 
        PlotStyle -> Hue[0], 
        DisplayFunction -> Identity]
lp2 = 
      ListPlot[eigvec[[2]], PlotJoined -> True, 
        PlotStyle -> Hue[1/3], 
        DisplayFunction -> Identity]
lp3 = 
      ListPlot[eigvec[[3]], PlotJoined -> True, 
        PlotStyle -> Hue[2/3], 
        DisplayFunction -> Identity]
Show[lp1, lp2, lp3, 
      DisplayFunction -> $DisplayFunction]

Now, let's make a Do loop that solves the problem for the potential V = (r^2) + ϵ (r^4) for ϵ = 0  …1  .0:

elist = {}; 
  

    vecplotlist = {};

    Do[
h   =  0.05; 
alku  =  
        0; 
loppu =  
        8; 
nn = (loppu - alku)/h + 
          1; 
Clear[x, r, v, i, j, 
        d]; 
v = 
        Table[r^2 + eps r^4, {r, alku, loppu, 
            h}]; 
d = -1/h^2; 
mat = 
        Table[Switch[i - j, -1, d, 0, v[[i]] -  2 d, 1, d, _, 
            0], {i, nn}, {j, nn}]; 
eigen = 
        Reverse[Eigenvalues[mat, -10]]; 
eigvec = 
        Reverse[Eigenvectors[mat, -10]]; 
elist = 
        Append[elist, 
          Table[{eps, eigen[[i]]}, {i, 1, 3}]]; 
lp1 = 
        ListPlot[eigvec[[1]], PlotJoined -> True, 
          PlotStyle -> Hue[0], 
          DisplayFunction -> Identity]; 
lp2 = 
        ListPlot[eigvec[[2]], PlotJoined -> True, 
          PlotStyle -> Hue[1/3], 
          DisplayFunction -> Identity]; 
lp3 = 
        ListPlot[eigvec[[3]], PlotJoined -> True, 
          PlotStyle -> Hue[2/3], 
          DisplayFunction -> Identity]; 
vecplotlist = 
        Append[vecplotlist, Show[lp1, lp2, lp3]], 
{eps, 0, 
        1.0, 0.1}]

Next, show all plots!

ListPlot[Flatten[elist, 1]];

nplot = 
        Length[vecplotlist];

    Table[
      Show[vecplotlist[[i]], DisplayFunction -> $DisplayFunction, 
        PlotRange -> {-0.25, 0.25}], {i, 1, nplot}]