Sciencemadness Discussion Board - xy” = (ax + b)y ?

xy” = (ax + b)y ?

blogfast25 - 8-4-2011 at 13:11

Is anyone here good with unusual differential equations? Then try this one for size!

A second order, non-linear DE:

xy” = (ax + b)y with a and b constants (Eq.1).

I’ve tried a few tricks, but no luck so far.

1. Reduce to 1st degree DE:

Substitute with y’ = zy

So y” = z’y + zy’ = z’y + z^2 and y” = y(z’ + z^2)

And x(z’ + z^2) = ax + b. This is 1st degree but second order in y DE and I can’t solve it. Shame, because with z = y’/y the second integration should be easy!

2. Assume a certain shape for y:

The real life problem that this DE represents states that for x = 0, y = 0 and for x = ∞, y = 0. This suggests a shape for y as follows:

y = f(x).e^(-kx)

If b = 0 then Eq.1 reduces to y” - ay = 0, a second order, linear, homogenous DE with constant coefficients and with a particular integral for k = √a, i.e. y = e^(-√a .x) for the homogeneous DE, or:

y = f(x).e^(-√a .x) seems a plausible shape if f(x) has a root at x = 0

It can be shown that if y = f(x).e^(-kx), then y” = g(x).e^(-√a .x) (Eq.2)

And g(x) = f”(x) - 2√a f’(x) + a f(x)

And thus using Eq.1, then x.g(x).e^(-√a .x) = (ax + b).f(x).e^(-√a .x)

Or x.g(x) = (ax + b).f(x)

So x.[ f”(x) - 2√a f’(x) + a f(x) ] = (ax + b) f(x) (Eq.3)

Trying f(x) = α x^2 + β x = x (α x + β

which has a root for x = 0 as required above:

f’(x) = 2 α x + β

f”(x) = 2α

Plug into Eq.3:

x.[2α - 2√a (2 α x + β

+ a.(α x^2 + β x)] = x (α x + β

(ax + b), drop x on both sides.

Left hand side = a α x^2 + (- 4 α √a + βa) x + (2α - 2√a β

Right hand side = a α x^2 + (α b + β a) x + β b

Comparison of coefficients:

For x^2: a α = a α

For x: - 4 α √a + β a = α b + β a or α = 0

And 2α - 2√a β = β b or with α = β = 0

Obviously f(x) = α x^2 + β x = x (α x + β

doesn’t fit.

But I still think some f(x) in the form of a polynomial with root x = 0 should found.

And this is as far as I got!

kmno4 - 11-4-2011 at 00:02

I would try this -
xy” = (ax + b)y can be writtren as homogenous equation y” - ([ax + b]/x)y = 0. Of course it is valid for x &#8800 0
Further steps - for example here:

Code:

http://mathworld.wolfram.com/EulerDifferentialEquation.html

woelen - 11-4-2011 at 00:45

What you have here is a so-called time-dependent homogeneous linear equation with (y y') being the state of the system and x being the time coordinate.

You can rewrite this into the form Y' = A(x)*Y, where Y is the state vector and A(x) a time-dependent coefficient matrix.

In general this kind of systems can only be solved analytically when the order of the system equals 1 or when A(x) is a constant matrix (i.e. does not depend on x). Your system has order 2 and is time-dependent. There is no general method for solving this kind of systems analytically and you'll have to fall back to either a numerical approximation of the solution (you need to know initial conditions then for y(0) and y'(0)) or you have to find by trial and error an analytic solution for this equation, but I doubt you will find one.

You might also have a more complicated problem when the fixation is not only at time 0, but you have values for y(x1) and y(x2) or y(x1) and y'(x2). In the latter case you have a boundary value problem which also needs to be approximated numerically.

Your problem actually is of the latter kind, you have y(0) = 0 and y(∞ ) = 0. By means of a coordinate transformation you can remove the seemingly singular nature of your problem and transform it into a boundary problem which does not involve infinity.

[Edited on 11-4-11 by woelen]

blogfast25 - 11-4-2011 at 05:36

kmno4:

The equation basically started off like that but re-writing it in that way doesn’t solve anything:

y” - (a + b/x)y = 0 (is of the general form y” + u(x).y’ + v(x).y = 0 but with u(x) = 0))

… is also a linear, second order homogeneous DE with non-constant coefficients. There is, as woelen states, no general mathematical solution for this type of equation, with the exception of a few ‘shapes’ like the Euler Cauchy DE, the Exact DE and probably a limited number of others that I’m not aware of.

Unless two particular solutions y1(x) and y2(x) are already known a priori, in which case the general integral can be found from y = c1.y1(x) + c2.y2(x) and the integration constants can then be determined from the boundary conditions. But that’s not our case either: I have trouble ‘guessing’ just one particular integral, never mind two!

Quote: Originally posted by woelen

Your problem actually is of the latter kind, you have y(0) = 0 and y(∞

= 0. By means of a coordinate transformation you can remove the seemingly singular nature of your problem and transform it into a boundary problem which does not involve infinity.

Switch to a coordinate system that removes the infinity? How would you do that?

woelen - 11-4-2011 at 07:09

I'll try to find a suitable transformation tonight.

In your attempt to solve this equation I see an error.

Asymptotic behavior for this equation (for x becoming very large) is

y'' = ay

This is in contradiction with your boundary condition y(∞ ) = 0 for any real a. For a > 0, the solution of this equation will diverge exponentially whatever the value of b. For a < 0, the solution of this equation will be oscillatory for very large x. But it might be that in such a case the oscillation is damped (very lightly) and that in the long run the amplitude goes to 0.

blogfast25 - 11-4-2011 at 08:06

woelen:

Yes, I see that now. Thanks. But y = x.(αx + β

.e^(-√a.x) doesn’t fit the DE anyway, as shown above…

Also, I’m sure that y is a symmetrical function with respect to x = 0: y(x) = y(-x) or at least |y(x)| = |y(-x)|, not sure which of the two. If true, it may be better to look for a solution with a Gaussian kernel in it: y = f(x). e^(-λ.x^2). And e^(-λ.x^2) descends much quicker too. Of course all this might be none other than a wild goose chase!

woelen - 11-4-2011 at 23:55

Asymptotically (x going towards +/- infinity) the solution can indeed be symmetric, but if this is the case, then the solution also is divergent and goes to infinity for both very negative x and very positive x.

Assume |x| is very large, then the equation simplifies to

y'' = ay

The solution of this equation is y = k1 * exp(-x *sqrt(a)) + k2 * exp(x * sqrt(a))

For positive real a, the solution is asymptotically symmetric if k1 == k2 and then the graph looks like a cosh() function for large |x|. In any case if you want y(∞ ) = 0 then k2 must be zero and then the asymptotic symmetry is lost and y(-∞ ) = ∞ .

So, if you want to find a solution with y(∞ ) = 0, then you have to look for a solution of the form

y = k1(x) * exp(-x *sqrt(a)), where the function k1(x) equals 0 at x = 0 and this function k1(x) approaches a non-zero constant for large |x| (here 'large' means large with respect to the absolute value of b).

From the properties of k1(x) you can also conclude that this is not a polynomial function, because every non-constant polynomial function goes to infinity when x goes to infinity.

Given the above analysis one can also take the cobnclusion that your guess with a Gaussian kernel will not lead to a solution. The nature of the solution simply is not of that form.

-----------------------------------------------------------------

There is one more question I have. Can you tell anything about the sign of a and b? If a and b have the same sign, then for x > 0 there is no singularity, otherwise there is a nasty thing at -a/b. If both are negative, then the asymptotic shape is not cosh()-like, but is oscillatory. If both are positive, then the asymptotic shape is cosh()-like or exponential-like.

[Edited on 12-4-11 by woelen]

blogfast25 - 12-4-2011 at 04:30

Woelen:

Thanks for your continued interest in this interesting little conundrum!

It’s probably useful at this point to reveal the actual ‘real world’ problem. It’s a little corker that my 16 year old daughter (now a student of science!) came up with (as only an innocent 16 year old can!): ‘What would a one dimensional hydrogen atom look like, dad?’ This came about after we’d been playing around with a Java applet for one-dimensional boxes, one-particle quantum systems, with various potential energy functions, their wave function and their energy levels. So I set out to set up the DE but it turned out dad had bitten off a little more than I could chew! Firstly, back to basics, to find the signs of a and b…

The wave function of a quantum system according to the time-independent form of the Schrödinger Equation fits:

Hψ = Eψ

With ψ the amplitude of the wave function, H the Hamiltonian Operator (HO) and E the total energy of the quantum system (Ekin + Epot).

For the one dimensional hydrogen atom the HO becomes:

- δ ψ(r)” - ρ ψ/r = E ψ

With:

• r the distance between the electron and the proton (assume the latter to be completely stationary)
• ψ (i.e. ψ(r)) the amplitude of the wave function as a function of r
• - ρ ψ/r the potential energy term resulting from the coulombic attraction between the electron and the proton
• δ and ρ two material constants that can be found in any decent physics book. Both are positive.
• E the total energy which in this case is always negative (E < 0), reaching 0 only for the unbound electron (r = ∞

Dividing everything by - δ and then multiplying by r gives:

r ψ” + (ρ/ δ

ψ = - (E/ δ

ψ r

or r ψ” = [- (E/ δ

r - (ρ/ δ

] ψ

Substituting: ψ = y, r = x:

xy” = (ax + b)y

with a = - (E/ δ

> 0
and b = - (ρ/ δ

< 0

Quote: Originally posted by woelen

So, if you want to find a solution with y(∞

= 0, then you have to look for a solution of the form

y = k1(x) * exp(-x *sqrt(a)), where the function k1(x) equals 0 at x = 0 and this function k1(x) approaches a non-zero constant for large |x| (here 'large' means large with respect to the absolute value of b).

From the properties of k1(x) you can also conclude that this is not a polynomial function, because every non-constant polynomial function goes to infinity when x goes to infinity.

Well, now I’m not so sure about symmetry anymore (or even whether it’s relevant) because x is a distance and can by its very definition only be x >= 0.

Also I’m not sure about x = 0, y = 0 anymore either.

Where I disagree is where you write: ”From the properties of k1(x) you can also conclude that this is not a polynomial function, because every non-constant polynomial function goes to infinity when x goes to infinity.”

Check out these functions for instance:

y = (2 - 2 √a . x) e^(- √a . x) and y = (27 - 54 √a . x + 18 a x^2) e^(- √a . x)

Both asymptote to y = 0 for x = ∞. The first one has one real root (x = 1/√a), the second two (x = 0.63/√a and x = 2.37/√a), so they ‘oscillate’ around y = 0, before settling for y = 0 for x = ∞. The polynomials used here are of the Laguerre type, like those that pop up in the solution for an actual (three dimensional) hydrogen atom. But the functions do not zero at x = 0...

I believe the solution can only be found by looking closer at the properties of the ‘real world’ system. For instance, if we go back to the ‘unperturbed’ (no potential energy term, or b = 0, free particle) wave function, y” - ay = 0, then as noted we get y = k1 e^(- √a . x)

And we know that a = - (E/ δ

and that for the bound electron (0 < r < ∞

the total energy E is quantised. So the sooner we can get a quantum number in there, the better. I believe the quantum number will be one of the integration constants, the second integration constant I believe will be 0.

We also know from the general properties of quantum systems that for the ground state E0, the wave function is entirely positive for 0 < r < ∞.

I hope this helps generate some ideas: I hate to disappoint my daughter!

EdiT: I've disabled smilies but the damn things keep popping up anyway, at least in my browser!

[Edited on 12-4-2011 by blogfast25]

watson.fawkes - 12-4-2011 at 05:23

Quote: Originally posted by blogfast25

And we know that a = - (E/ δ

and that for the bound electron (0 < r < ∞

the total energy E is quantised. So the sooner we can get a quantum number in there, the better

We know it's quantized in three dimensions. That doesn't mean it's quantized in one dimension. I doubt it is. Quantization for dimensions n >= 2 arises from phase matching on the sphere of dimension n-1. Since S0 is just a pair of points (+1, -1), it doesn't generate any such continuity requirement. Also notice that your Hamiltonian doesn't have an angular momentum term; this is related.

blogfast25 - 12-4-2011 at 07:26

Hmmm… good point, Watson. But the quantum systems for one-dimensional finite wells, infinite wells, harmonic oscillators and various potential energy functions are all quantised in 1 dimension. I’ve seen one for a coulomb but if I recall well that was with infinite well boundaries (Epot = - δ/r for -L < x < +L, Epot = ∞ outside the ‘box’), it too was quantised.

These have no angular momentum terms either. For the actual hydrogen atom the radial component R(r) is solved separately from the angular components and only the radial component is quantised.

[Edited on 12-4-2011 by blogfast25]

watson.fawkes - 12-4-2011 at 19:20

For ease of notation, I'll write the equation with a2 instead of a; this avoids all the square roots. With this change, the differential operator x D2 - a2 x + b has a zero (solution to the differential equation) when b = 2a, and the solution is f = x e-ax. This is the ground state solution. The first excited solution happens when b = 4a, and it is f = ( 1 - ax ) x e-ax. (Assuming I've done all my arithmetic right.)

I didn't get these solution by magic. I got out my old QM textbook and read up how the principal quantum number gets quantized. The differential equation in question is essentially that in the general case, but with l = 0. The argument is, to my adult mind, extraordinarily unsatisfying mathematically, but boils down to a power series expansion of an unknown function, deriving a recurrence relationship on the coefficients, and showing that the solution diverges unless the polynomial is of finite degree. The degree of the polynomial is the principal quantum number, a fact that I had forgotten. That means there's a solution of the form F(x) x e-ax, where F is a polynomial. Substituting that solution back into the original equation yields a polynomial with coefficients that are expressions in a, b, and c that must be identically zero to solve the equation. A polynomial is only identically zero if all of its coefficients are. That leads to both the relation between a and b and the coefficients of F. For deg F = 0, there's just the relation; the single coefficient can be taken as 1, since the differential equation is homogeneous.

Edit: Perhaps-whoops! Math error. Maybe. Verifying now.

[Edited on 14-4-2011 by watson.fawkes]

blogfast25 - 13-4-2011 at 08:21

Watson:

Firstly I agree with your comments regarding how weird the solution really is, with the quantum number popping up as the degree of the polynomial. But: t’is also thus for the solution of the actual hydrogen atom.

Regards the proposed solution, in my evaluation it’s close but no (full) cigar. In short it works for the ground state, but not for the first excitation. Let me explain.

The original DE was:

xy” = (ax + b)y

with a = - (E/δ ) > 0
and b = - (ρ/δ ) < 0

I’ve rewritten this once again as:

xy”= (α2x - β )y

with α2 = - (E/δ )
and β = - b = ρ/δ

This is of course the same equation as above

You proposed β = 2 α (in my notation) for the ground state.

Let’s look at the ground state energy:

α = β/2

Square both sides:

α2 = β2/4

which works out as E1 = - ρ3/(4 δ2

which is indeed ≠ 0 and negative as it should be.

The DE can then be re-written as xy”= (α2x - 2α )y

The proposed solution to the DE is y = xe- αx and it does indeed fit the equation xy”= (α2 - 2α )y. I’ve plugged y = xe- αx into it and it works!

Let’s now look at the first excitation:

You proposed (in my notation) α = β/4, so
E2 = - ρ3/(16 δ2

The factors 4 and 16 would suggest quantisation as (2n)2 with n the quantum number, so

En = - ρ3/((2n)2 δ2

The DE for n = 2 would then be:

xy”= (α2x - 4α )y

And you proposed y2= x(1 - αx)e- αx but that doesn’t fit the DE.

Twice differentiating, multiplying by x and dividing by e- αx I get for the left hand side:

Edit: a small correction was carried out here.

-4αx + 4α2x2 + α2x2 - α2x3

And for the right hand side:

α2x2 - α3x3 - 4αx + 4α2x2

So you end up with:

4α2x2 - α2x3 = -α3x3 + 4α2x

Two possibilities: either there’s an error in my differentiations (I will rigorously check again tonight) or the polynomial x(1 - αx) is incorrect. I haven’t tried any others yet. Note that I’d already tried a polynomial of that kind above and couldn’t get it to work either.

In the case the polynomial is incorrect it’ll be a case of ”but boils down to a power series expansion of an unknown function, deriving a recurrence relationship on the coefficients, and showing that the solution diverges unless the polynomial is of finite degree”. That should be fun [cough!]

[Edited on 13-4-2011 by blogfast25]

[Edited on 13-4-2011 by blogfast25]

watson.fawkes - 13-4-2011 at 09:44

Use the expression (1 + cx) x e-ax as a candidate solution. You'll get some simple expressions for the coefficients of a polynomial. All of these coefficients must be zero.

I may have made a clerical error doing the algebra by hand.

smaerd - 13-4-2011 at 10:02

I have nothing to add here except wow, differential equations looks a lot harder then calculus 2!

blogfast25 - 13-4-2011 at 12:10

Smaerd:

This is a notoriously hard one. Well above my pay grade!

Quote: Originally posted by watson.fawkes

Use the expression (1 + cx) x e-ax as a candidate solution. You'll get some simple expressions for the coefficients of a polynomial. All of these coefficients must be zero.

I may have made a clerical error doing the algebra by hand.

Regards the clerical error, I think I spotted one high up and corrected it without passing it down.

I also wonder if we’re overlooking something: if you look at the exponential part of a radial wave function for a monoelectric atom you’ll see that the exponent of e- αx is in fact not constant but varies with the quantum number. See here for example (look at the exponential function: its exponent depends on n):

http://panda.unm.edu/Courses/Finley/P262/Hydrogen/WaveFcns.h...

Well, assuming the energy quantisation is indeed En = - ¼ (ρ3/δ2

(1/n)2, then it can be shown easily with α2 = - (E/δ ) that α = ½ (ρ/δ )3/2(1/n)

So I’m wondering if looking for a function:

yn = f(x,n)e^[- ½ (ρ/δ )3/2(1/n) x]

… would shed any light on the matter, with f(x,n) a polynomial. These Laguerre polynomials keep calling on me…

It’d be nice to at least solve the first excitation because I think it would prove that En = - ¼ (ρ3/δ2

(1/n)2, is correct.

[Edited on 13-4-2011 by blogfast25]

watson.fawkes - 13-4-2011 at 19:04

So I got out ye olde computer algebra. Here's a brute force expression that illustrates how to find the solution.<pre>In:= x D[(1 + c x) x Exp[-a x], x, x] - a^2 x (1 + c x) x Exp[-a x] + b (1 + c x) x Exp[-a x]
Out:= E^(-a x) x (-2 a + b + 2 c + (-4 a + b) c x)</pre>From the coefficient of x, you get the relation b = 4a. From the constant coefficient, you get c = -a.

Looking at your computation of the right hand side, the last term should be a3 x3, not a2 x3. (Verified with software.) Then the two sides match.

Just for fun, I can brute-force the second excited state.<pre>In:= x D[(1 + c x + d x^2) x Exp[-a x], x, x] - a^2 x (1 + c x + d x^2) x Exp[-a x] + b (1 + c x + d x^2) x Exp[-a x]
Out:= E^(-a x) x (-2 a + b + 2 c + (-4 a c + b c + 6 d) x + (-6 a + b) d x^2)</pre>OK. That's the messy work. What's the polynomial and the value of b?

woelen - 13-4-2011 at 23:24

Good to see that there is a solution now to the physical problem, but it is not a solution to the original problem for general a and b. Here we have a solution for specific combinations of a and b, derived from quantum states.

I assume the solution for general a and b is not relevant anymore.

The fact that the quantum number pops up inthe degree of a polynomial is not that strange to my opinion. There are so many physical processes,, where the order of the process appears as a polynomial of that order.

Also look at the orbitals of electrons around atoms and look at their shapes. They start of as simple spheres (e.g. depicting the 90% chance of presence volume) and become progressively more complicated. The description of these figures also requires more and more complicated formulas (polynomials) for increasing order.

blogfast25 - 14-4-2011 at 05:51

Watson:

Thanks a bunch! As it so happens, yesterday evening I plugged in y = x(1- cx ) e- αx into xy” = (α2x - 4α )y (once again!)

And obtained… c = α !

So, my bad and we’re now in full agreement: I also found the original error which was banal and typical of fairly longish manual algebra (what digital gizmo are you using for this brute force application??) Once you make this type of error you tend to repeat it like a demented parrot, it becomes ‘taught behaviour’. Very annoying (but ‘never again’!)

So y = x(1- αx ) e- αx is a solution for β = 4α and presumably other solutions exist for β = 6α, β = 8α etc.

Re. quantisation:

En = - ρ3/((2n)2 δ2

And from above α2 = - (E/δ ), E = - α2 δ

Thus: αn = ½ (ρ/δ )3/2 (1/n)

Set α1 = ½ (ρ/δ )3/2

Then α2 = α1/2

α3 = α1/3

...

αn = α1/n

So y2 = x (1 - ½ α1 x ) e^(- ½ α1 x)

Which makes sense I guess (w/o having plotted anything yet), because at higher energy the electron will venture further away from the ‘nucleus’.

I'll have a look at your data for the second excitation a bit later on, I think β = 6α.

Can we now perhaps try and envisage a general solution?

Say yn = x F(x, n, αn

e^(- αnx)

With F = 1 + c1x + c2x2 + … + cn-1xn-1

Or F = [1+ (i from 1 to n-1) ∑ cixi]

And then there’s the ‘minor’ problem of normalisation: will yn have to be normalised so that (from 0 to ∞ ) A2 ∫ y2dx = 1 (A is the normalisation constant)? I’m guessing ‘yes’ and will have a go at y1 later on...

Woelen:

Quote: Originally posted by woelen

I assume the solution for general a and b is not relevant anymore.

Also look at the orbitals of electrons around atoms and look at their shapes. They start of as simple spheres (e.g. depicting the 90% chance of presence volume) and become progressively more complicated. The description of these figures also requires more and more complicated formulas (polynomials) for increasing order.

No, the general solution remains relevant. But the clue is in the last part of your post: basically to derive a general solution we’d have to apply similar techniques as applied to the radial part of the hydrogen wave function. That math is somewhat impenetrable for my level but I’ve found some good texts about and will ‘have a go’, time permitting!

[Edited on 14-4-2011 by blogfast25]

[Edited on 14-4-2011 by blogfast25]

watson.fawkes - 14-4-2011 at 06:18

Quote: Originally posted by blogfast25

Can we now perhaps try and envisage a general solution?

You've got the basic idea, which is to write a candidate solution in a power series and then plug it into the differential equation. Each coefficient of this series must be zero, so it gives a series of relationships between coefficients. This equation is only second order in derivatives, which means the recurrence relationship is at most second order. For example, the defining recurrence of the Fibonacci series is second order, as it depends on two previous terms. As it turns out, the recurrence derived from this equation is single order.

The right expression is x F(x) e-ax, and write F(x) = &Sigma cixi, taking c0=1.

blogfast25 - 14-4-2011 at 08:26

Watson:

Agreed.

And it’s official, the second excitation is y3 = x(1 - 2αx + 2/3.α2x2

e- αx

With β = 6α and α = α1/3

blogfast25 - 15-4-2011 at 09:09

Setting up the polynomial using F = 1 + ∑cixi (from i = 1 to n - 1) was surprisingly easy. First y was transformed to y = Ge- αx
with G = x + ∑cixi+1.

y” can be written as e- αx [G” - 2α G’ + α2 G]

The x and the e- αx in the equation xy” = (α2x - 2n α )xF e- αx drop out and the equation becomes:

G” - 2α G’ + α2 G = (α2x - 2n α )F

After deriving G twice and a bit of rearranging, the full polynomial is then obtained. All its coefficients need to be zero to satisfy the DE. The recursing is real easy because it turns out that ci is always a simple function of ci-1:

ci = 2/(i(i+1)) . α (i-n) ci-1 and c0 = 1

For i=1, c1 = α (1-n) which fits for n=1 (c1 = 0) and n=2 (c1 = - α ) and for n=3 (c1 = - 2α ), all as found above.

For i=2, c2 = 1/3 α2(2-n)(1-n)

Which for n=3 gives c2 = 2/3 α2, as found above.

c3 = 1/18 . α3(3-n)(2-n)(1-n)

Etc etc.

Note also that in all cases α = α1/n

I could stop there but a general expression for yn = U(x, n) will be needed because the functions need normalising. For y1 the normalisation constant is 2 α13/2. Without normalising it will be impossible to plot comparable graphs for n=1, n=2… n

[Edited on 15-4-2011 by blogfast25]

watson.fawkes - 15-4-2011 at 11:27

Quote: Originally posted by blogfast25

ci = 2/(i(i+1)) . α (i-n) ci-1 and c0 = 1

It's the i - n term that causes the series to terminate at some point. It has to be zero, and then every subsequent coefficient is also zero. You can substitute b later in the process, and thereby derive permissible ratios for b /a (they're all rational).

So I take it that "a power series expansion of an unknown function, deriving a recurrence relationship on the coefficients" doesn't seem quite so daunting now? It's much easier after you see the low-order specific cases first.

blogfast25 - 16-4-2011 at 04:23

Yes, I can see how even for higher degree DEs of that type, say y(n)' = f(x)y with f(x) a simple function) this method with anticipated solutions y = (polynomial).e- kx can work because the degree of the polynomial isn’t affected by the number of derivations. I can even see how the method could be generalised for some specific cases. Very interesting stuff!

It should also work for y(n)’+ y(n-1)’ + y(n-2)’ + … + y" + y’ forms, I’m guessing…

But what I still don’t really understand is how did you ‘guess’ or work out that for our case there were solutions for β = 2 n α? That α was quantised is a given because it’s basically the total energy term. But the exact form of the quantisation?

[Edited on 16-4-2011 by blogfast25]

watson.fawkes - 16-4-2011 at 07:34

Quote: Originally posted by blogfast25

But what I still don’t really understand is how did you ‘guess’ or work out that for our case there were solutions for β = 2 n α? That α was quantised is a given because it’s basically the total energy term. But the exact form of the quantisation?

Well, the short answer is that I didn't. I looked up the solution in my old QM textbook. I fed it back in bite-size pieces.

Choosing the general form to examine was done in the following way. First examine the solutions in the limit when x --> ∞. The solution to this is the term e-ax. Then approximate the solutions in the neighborhood of x = 0. The solution here is x. Dividing out by these factors, the remaining factor in a solution (a) has any value at x = 0, and (b) is sub-exponential at infinity. What this does is make a more sane recurrence on the coefficients. Luckily, it's sub-exponential by being polynomially-bounded (being a polynomial itself). It could also be sub-exponential in other ways, such as e√x, in which you are unlucky and have to do more work.

The final answer to "how", though, is only justified after the fact. As with all differential equations, there's no general algorithm.

blogfast25 - 16-4-2011 at 13:06

I see.

The general formula for ci appears to be:

ci = 2i/i!(i+1)!) αi Π ((i-n)(i-1-n)(i-2-n)…(2-n)(1-n))

with Π from i to 1 and α= αi/n

Inserting this into G promises lots of [ahem!] fun for the normalisation, which requires the determined integral (0 to ∞ ) ∫ G2e^(- 2 α1x) dx to be calculated. Well, for the simpler ∫ xm ex dx the integral is recursive: ∫ xm ex dx = xm ex - m ∫ xm-1 ex dx + C, if I recall well.

[Edited on 16-4-2011 by blogfast25]

blogfast25 - 18-4-2011 at 08:20

I decided to bypass full normalisation and replace it with numerical normalisation (integrating the functions y2 by means of a spread sheet, used also for the plots).

Below are the normalised wave functions obtained that way, with α1 = ½ (2π e2m/hε0

3/2 (e = charge of the electron, m = mass of the electron, h = Planck’s constant, ε0 = vacuum permittivity) - ground state and 1st and 2nd excitations:

The electron densities:

Note how the second and third 'lobes' are always higher than the first one.

[Edited on 18-4-2011 by blogfast25]