Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am in the process of trying to use Mathematica to find the coefficients of transmission and reflection for the wavefunction of a beam of particles tunnelling across any number of identical potential barriers.

I have got to the stage where I am able to generate the equations for n numbers of barriers, but I am not able to equate the appropriate equations to each other to solve them for the coefficients. The module I have written to generate the relavant equations is:

GeneratePhiEqs[N_] :=
    Module[{inceqs, diffinceqs, bareqs, diffbareqs, finaleq, difffinaleq,
    IncidentPhi, BarrierPhi, FinalTransmittedPhi, phieqs}, 
    IncidentPhi[x_, n_] := T[n]*E^(I*k*x) + R[n]*E^(-I*k*x);
    BarrierPhi[x_, n_] := A[n]*E^(r*x) + B[n]*E^(-r*x);
    FinalTransmittedPhi[x_, n_] := T[n]*E^(I*k*x); T[0] = 1;

    inceqs = Table[Phi[j] = {IncidentPhi[x, j]}, {j, 0, N}];
    diffinceqs = 
    Table[Phi[N + j] = {D[IncidentPhi[x, j], x]}, {j, 0, N - 1}];
    bareqs = Table[Phi[2 N + j - 1] = {BarrierPhi[x, j]}, {j, N}];
    diffbareqs = 
    Table[Phi[3 N + j - 1] = {D[BarrierPhi[x, j], x]}, {j, N}];
    finaleq = Phi[4 N] = FinalTransmittedPhi[ x, N];
    difffinaleq = Phi[4 N + 1] = D[FinalTransmittedPhi[x, N], x];
    phieqs = {inceqs, diffinceqs, bareqs, diffbareqs, finaleq, 
    difffinaleq}; phieqs]

The module I have written (or what it looks like so far) to equate the relevant equations is:

SolvePhiEqs[number_] := 
    Module[{eqs, simeqs}, eqs = GeneratePhiEqs[number]; 
    simeqs = Table[{Phi[j] /. phieqs == Phi[2*number + j] /. 
    phieqs}, {j, 2*number - 1}]; simeqs] 

Basically, the equations that I generated with the first module will not substitute in correctly, they just remain isolated from each other surrounded by curly brackets.

Any suggestions would be much appreciated!

Thanks.

share|improve this question
    
I get a collection of stuff by running GeneratePhiEqs[1]. But these are not equations, yet, they have no equality. Can you explain what terms are to be equated with what? –  bill s Jul 10 '13 at 13:38
    
Can you write what explicitly A[n], B[n], R[n], T[n] are? Also, the usage Table[Phi[j] = {IncidentPhi[x, j]}, {j, 0, N}] looks so weird to me. I would rather go with something much simpler like Phi[j_]:={IncidentPhi[x, j]}. –  Ali Jul 10 '13 at 14:03
    
@bills GeneratePhiEqs effectively gives me a set of functions and their derivatives, some of which I know are equal to each other at a certain value of x. I need to equate particular functions to each other and substitute in a value for x in order to solve for the coefficients. The SolvePhiEqs module should take care of equating them to each other and solving them. –  CJO Jul 10 '13 at 14:23
    
In answer to your question @Ali the A[n], B[n] are coefficients. For example, if I run GeneratePhiEqs[1], I get a certain number of equations, and it numbers the coefficients for me. I am trying to solve for the coefficients. –  CJO Jul 10 '13 at 14:24
    
I get a function like $e^{ikx}$ at one barrier, but then I get the same thing at another barrier. I want to differentiate between them whilst still being able to see that they are related. –  CJO Jul 10 '13 at 14:29
show 1 more comment

1 Answer

up vote 3 down vote accepted

Possibly this?

Clear[Phi];
GeneratePhiEqs[N_] := (* <- Don't use N *)
 Module[{inceqs, diffinceqs, bareqs, diffbareqs, finaleq, difffinaleq,
    IncidentPhi, BarrierPhi, FinalTransmittedPhi}, 
  IncidentPhi[x_, n_] := T[n]*E^(I*k*x) + R[n]*E^(-I*k*x);
  BarrierPhi[x_, n_] := A[n]*E^(r*x) + B[n]*E^(-r*x);
  FinalTransmittedPhi[x_, n_] := T[n]*E^(I*k*x); T[0] = 1;
  inceqs      = Table[Phi[j] -> IncidentPhi[x, j], {j, 0, N}];
  diffinceqs  = Table[Phi[N + j] -> D[IncidentPhi[x, j], x], {j, 0, N - 1}];
  bareqs      = Table[Phi[2 N + j - 1] -> BarrierPhi[x, j], {j, N}];
  diffbareqs  = Table[Phi[3 N + j - 1] -> D[BarrierPhi[x, j], x], {j, N}];
  finaleq     = Phi[4 N] -> FinalTransmittedPhi[x, N];
  difffinaleq = Phi[4 N + 1] -> D[FinalTransmittedPhi[x, N], x];
  Flatten @ {inceqs, diffinceqs, bareqs, diffbareqs, finaleq, difffinaleq}]

SolvePhiEqs[number_] := 
 Module[{eqs},
  eqs = GeneratePhiEqs[number];
  Table[Phi[j] == Phi[2*number + j] /. eqs, {j, 2*number - 1}]
 ]

If it's right, then you should

Example:

SolvePhiEqs[2]
(* {E^(-I k x) R[1] + E^(I k x) T[1] == E^(r x) A[2] + E^(-r x) B[2], 
    E^(-I k x) R[2] + E^(I k x) T[2] == E^(r x) r A[1] - E^(-r x) r B[1],
    -I E^(-I k x) k R[1] + I E^(I k x) k T[1] == E^(r x) r A[2] - E^(-r x) r B[2]} *)
share|improve this answer
    
It works - Thank you so much! Please excuse my cobbled together syntax, I'm a bit of a Mathematica beginner. –  CJO Jul 10 '13 at 16:01
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.