Mathematica help (plotting solutions of a nonlinear equation)

[maverick280857]

To the moderator: Please shift this to an appropriate forum, if necessary.

Hi,

I am solving a QM problem which requires using Mathematica to solve two equations,

[tex]J_{\alpha}'(2\sqrt{\lambda}) = 0[/tex]

(derivative of Bessel function)

and

[tex]J_{\alpha}(2\sqrt{\lambda}) = 0[/tex]

(not simultaneously) where [itex]\alpha[/itex] is the unknown and [itex]\lambda[/itex] is known. I want to plot the dependence of [itex]\alpha^{2}[/itex] on [itex]\lambda[/itex].

First, I take [itex]x = 2\sqrt{\lambda}[/itex] and vary [itex]x[/itex] from 1 to 10 (so [itex]\lambda \in [0.25, 25][/itex]). Now, I use

FindRoot[0.5*(BesselJ(-1+alpha,x) - BesselJ(1+alpha,x)), {alpha,0}]

for different values of [itex]x[/itex] which are being keyed in by hand while typing this command. To automate it, I set up an array:

Array[SOL, 1000]

and then use a for loop

For[c=1;x=1,c<=1000;x<=100,c = c + 1; x=x+0.01, SOL[c] = {x,FindRoot[(BesselJ[-1+alpha,x] + BesselJ[1+alpha,x])/2, {alpha,0}]}]

so SOL[10] is of the form:

{1.09, {alpha -> 0.448962}}

I have the following questions:

1. The second element of the ordered pair above is of the form {alpha --> value}. How do I extract alpha from this? I want to store it in (x, alpha) form.

2. How do I plot some function of alpha^2 versus some function of x, from this data?

3. Is there an easier way to do all this?

I've been trying to figure this out for some time but I haven't had much success.

Thanks in advance.

Cheers
Vivek

 

[vanesch]

To the moderator: Please shift this to an appropriate forum, if necessary.

Hi,

I am solving a QM problem which requires using Mathematica to solve two equations,

[tex]J_{\alpha}'(2\sqrt{\lambda}) = 0[/tex]

(derivative of Bessel function)

and

[tex]J_{\alpha}(2\sqrt{\lambda}) = 0[/tex]

(not simultaneously) where [itex]\alpha[/itex] is the unknown and [itex]\lambda[/itex] is known. I want to plot the dependence of [itex]\alpha^{2}[/itex] on [itex]\lambda[/itex].

First, I take [itex]x = 2\sqrt{\lambda}[/itex] and vary [itex]x[/itex] from 1 to 10 (so [itex]\lambda \in [0.25, 25][/itex]). Now, I use

FindRoot[0.5*(BesselJ(-1+alpha,x) - BesselJ(1+alpha,x)), {alpha,0}]

for different values of [itex]x[/itex] which are being keyed in by hand while typing this command. To automate it, I set up an array:

Array[SOL, 1000]

and then use a for loop

For[c=1;x=1,c<=1000;x<=100,c = c + 1; x=x+0.01, SOL[c] = {x,FindRoot[(BesselJ[-1+alpha,x] + BesselJ[1+alpha,x])/2, {alpha,0}]}]

so SOL[10] is of the form:

{1.09, {alpha -> 0.448962}}

I have the following questions:

1. The second element of the ordered pair above is of the form {alpha --> value}. How do I extract alpha from this? I want to store it in (x, alpha) form.

If you have, say: u = {1.09, {alpha -> 0.448962}}
then you can obtain the alpha value by doing:
alpha /. u[[2]]

(the u[[2]] extracts the second element, namely {alpha -> 0.44...} and alpha /. applies this as a replacement rule on the expression "alpha".

2. How do I plot some function of alpha^2 versus some function of x, from this data?

The simplest thing is probably to construct first a list of the kind:

uu = Map[({#[[1]], alpha/. #[[2]] })&, SOL ]

or something similar (didn't check it on mathematica, hope it works more or less).

After that, you do a ListPlot of the array.

3. Is there an easier way to do all this?

It is an easy way!

 

[maverick280857]

Thanks vanesch :-)

If you have, say: u = {1.09, {alpha -> 0.448962}}
then you can obtain the alpha value by doing:
alpha /. u[[2]]

This works fine for u, but if I have an element of SOL, say SOL[2], then using the same idea:

alpha /. SOL[[2]]

or

alpha /. SOL[2]

both do not work.

 

[alphysicist]

Hi maverick280857,

If I'm understanding your question corrrectly, then if

SOL[[10]]= {1.09, {alpha -> 0.448962}}

you could use:

alpha /. SOL[[10,2]]

to extract alpha.



Or, inside your For loop, instead of:

... SOL[[c]]= { x, FindRoot [...] } ...

I think you could get around the issue by getting the alpha value directly:

... SOL[[c]]= { x, alpha/. FindRoot [...] } ...

 

[maverick280857]

Thanks alphysicist.

Hi maverick280857,

If I'm understanding your question corrrectly, then if

SOL[[10]]= {1.09, {alpha -> 0.448962}}

you could use:

alpha /. SOL[[10,2]]

to extract alpha.

This did not work either. And also, its SOL[10] which holds the value {1.09, {alpha -> 0.448962}} and not SOL[[10]].

Also, using SOL[[10]] gives the following message:

Part::partd: Part specification SOL〚10〛 is longer than depth of object.

 

[alphysicist]

Oh, I see what you have. Then you can just use:

SOL[10][[2]]



What about the idea of using

... SOL[c]= { x, alpha/. FindRoot [...] } ...

inside the FOR loop so that SOL[10] is what you want it to be in the first place?




By the way, I don't think that your commands are doing what you intend. The command

Array[SOL,1000]

is not setting up a 1000 element array named SOL. It is displaying a 1000 element array whose elements are the variable SOL, but it does not do anything with that array; it looks like it's being lost. You could run the For loop without running this line at all.

If you want to create a 1000 element array named SOL with all zero elements you might try:

SOL=Table[ 0, {i,1,1000}]

or you could use a DO or FOR loop. Then the tenth element of the array SOL is accessed by SOL[[10]].

In your For loop, when you use SOL[c] (with one bracket on each side) you are not accessing elements of an array but assigning values to a function SOL. (That's what the error meant in your last post: SOL is not an array.) Once you create the matrix SOL, you can then use SOL[[c]] inside your For loop to assign values to the matrix, and you can use:

...SOL[[c]]={x, alpha/. FindRoot[...]} ...

in the FOR loop to keep from getting the 'arrowed' solution.



Of course there are plenty of other ways, too. Mathematica has lots of options.

 

[maverick280857]

It worked this time. Thanks :-)

Now I have SOL of the form

[itex][x_{i}, y_{i}][/itex]

and I want to plot these points for [itex]1 \leq i \leq i_{max}[/itex]. How can I do this? Plot and ListPlot don't work.

 

[vanesch]

It worked this time. Thanks :-)

Now I have SOL of the form

[itex][x_{i}, y_{i}][/itex]

and I want to plot these points for [itex]1 \leq i \leq i_{max}[/itex]. How can I do this? Plot and ListPlot don't work.

You're making your life difficult by having put the value in SOL, and not in SOL[], because the first form is a kind of function definition for a single value, while the second is an element of a list. ListPlot can only plot lists, so you could maybe first create a dummy list
SOL = Range[1,1000] and then use SOL[[c]] = ...

 

[alphysicist]

maverick280857,

That's why I was referring to your method in my last post. SOL in your calculation is not an array, and so you are cutting yourself off from the list operations that mathematica has.

That's why ListPlot doesn't work: because SOL is not an array. You've set it up as a function. However, it's a function that you have only defined at specific points, and inbetween those points it's undefined. That's why Plot isn't working.


I would suggest again that you create SOL as an actual array (or list I guess they call it). I prefer using the Table command to create straightforward lists, so something like

SOL= Table[ {x, alpha/. FindRoot[(BesselJ[-1+alpha,x] - BesselJ[1+alpha,x])/2, {alpha,0}]}, {x,1,100,0.01}]

as the first command; this sets up your list so that x varies from 1 to 100, in increments of 0.01. You can then access the 10th element using SOL[[10]] which gives {1.09,0.448962}. If you want to plot it, you can use

ListPlot[SOL]


If you would rather use Plot instead of ListPlot, you can create an interpolating function to find the points between the given points and then plot it with the following two commands:

SOLFUNC=Interpolation[SOL]

Plot[SOLFUNC[x],{x,1,100}]