NDsolve problem
in [Mathematica]
|
From: ntina.86 on 2 Jul 2010 07:27 Hello, I haven't been using Mathematica for some years and now I need to ask you for help with NDSolve. I'm using version 6.0. I want to solve the radial Schrodinger equation with a Yukawa potential, and I started trying to solve this: a = 500 b = 10; s = NDSolve[{y''[x] + (2/x) y'[x] + (1 + (2*a/x)*E^(-b*x)) (y[x]) == 0, y [0.000001] == 1, y'[0.000001] == -a}, y, {x, 0.000001, 60},PrecisionGoal -> \ [Infinity], MaxSteps -> Infinity] This gives me results that might be correct, but the problem comes when I try to plot: Plot[Evaluate[(x*y[x])^2 + (((x - 0.5 Pi)*y[x - 0.5 Pi]))^2 /.s[[1]]], {x, 10, 60}, PlotRange -> All] It should be a constant (i.e. the sum of the square of a cosine and of a sine, which is the solution of the Schrod.) at least for x>10, instead I get a badly oscillating function. Is it a problem of the method of integration or I'm doing something wrong? I tried to change it, using ExplicitRK or the SymplecticPartitionedRK, but the latter gives me these errors: --------------------------------------------------------------------------------------------------------- NDSolve::nlnum1: The function value {-21065.3 (1.23594*10^7+82.8327 \ TemporaryVariable[2][1.1462*10^-6])} is not a list of numbers with \ dimensions {1} when the arguments are {1.1462*10^-6,0.494377}. InterpolatingFunction::dmval: Input value {28.4292} lies outside the \ range of data in the interpolating function. Extrapolation will be \ used. >> InterpolatingFunction::dmval: Input value {30} lies outside the range \ of data in the interpolating function. Extrapolation will be used. >> NDSolve::nlnum1: The function value {-21193.5 (1.23496*10^7+82.3232 \ TemporaryVariable[2][1.14631*10^-6])} is not a list of numbers with \ dimensions {1} when the arguments are {1.14631*10^-6,0.493983}. InterpolatingFunction::dmval: Input value {28.4292} lies outside the \ range of data in the interpolating function. Extrapolation will be \ used. >> General::stop: Further output of InterpolatingFunction::dmval will be \ suppressed during this calculation. >> NDSolve::nlnum1: The function value {-21321.9 (1.23397*10^7+81.8195 \ TemporaryVariable[2][1.14643*10^-6])} is not a list of numbers with \ dimensions {1} when the arguments are {1.14643*10^-6,0.493589}. General::stop: Further output of NDSolve::nlnum1 will be suppressed \ during this calculation. >> ------------------------------------------------------------------------------------------------------- Thank you in advance. Regards, Valentina
From: Peter Pein on 3 Jul 2010 08:19 Am Fri, 2 Jul 2010 11:27:28 +0000 (UTC) schrieb "ntina.86(a)libero.it" <ntina.86(a)libero.it>: > Hello, > I haven't been using Mathematica for some years and now I need to ask > you for help with NDSolve. > I'm using version 6.0. > > I want to solve the radial Schrodinger equation with a Yukawa > potential, and I started trying to solve this: > > a = 500 > b = 10; > s = NDSolve[{y''[x] + (2/x) y'[x] + (1 + (2*a/x)*E^(-b*x)) (y[x]) == > 0, y [0.000001] == 1, y'[0.000001] == -a}, y, {x, 0.000001, > 60},PrecisionGoal -> \ [Infinity], MaxSteps -> Infinity] > > This gives me results that might be correct, but the problem comes > when I try to plot: > > Plot[Evaluate[(x*y[x])^2 + (((x - 0.5 Pi)*y[x - 0.5 > Pi]))^2 /.s[[1]]], {x, 10, 60}, PlotRange -> All] > > It should be a constant (i.e. the sum of the square of a cosine and > of a sine, which is the solution of the Schrod.) at least for x>10, > instead I get a badly oscillating function. > > Is it a problem of the method of integration or I'm doing something > wrong? I tried to change it, using ExplicitRK or the > SymplecticPartitionedRK, but the latter gives me these errors: > > > --------------------------------------------------------------------------------------------------------- > NDSolve::nlnum1: The function value {-21065.3 (1.23594*10^7+82.8327 \ > TemporaryVariable[2][1.1462*10^-6])} is not a list of numbers with \ > dimensions {1} when the arguments are {1.1462*10^-6,0.494377}. > > InterpolatingFunction::dmval: Input value {28.4292} lies outside the \ > range of data in the interpolating function. Extrapolation will be \ > used. >> > > InterpolatingFunction::dmval: Input value {30} lies outside the range > \ of data in the interpolating function. Extrapolation will be used. > >> > > NDSolve::nlnum1: The function value {-21193.5 (1.23496*10^7+82.3232 \ > TemporaryVariable[2][1.14631*10^-6])} is not a list of numbers with \ > dimensions {1} when the arguments are {1.14631*10^-6,0.493983}. > > InterpolatingFunction::dmval: Input value {28.4292} lies outside the \ > range of data in the interpolating function. Extrapolation will be \ > used. >> > > General::stop: Further output of InterpolatingFunction::dmval will be > \ suppressed during this calculation. >> > > NDSolve::nlnum1: The function value {-21321.9 (1.23397*10^7+81.8195 \ > TemporaryVariable[2][1.14643*10^-6])} is not a list of numbers with \ > dimensions {1} when the arguments are {1.14643*10^-6,0.493589}. > > General::stop: Further output of NDSolve::nlnum1 will be suppressed \ > during this calculation. >> > > > ------------------------------------------------------------------------------------------------------- > > > Thank you in advance. > Regards, > Valentina > > > Hi Valentina, it's strange. With Mathematica 7 but without the Options for NIntegrate, I get a plot-function which is (almost) constant a little bit below 1/1000: In[1]:= a=500; b=10; \[CurlyEpsilon]=10^-6; s=NDSolve[{y''[x]+(2/x) y'[x]+(1+(2*a/x)*E^(-b*x)) (y[x])==0,y[\[CurlyEpsilon]]==1,y'[\[CurlyEpsilon]]==-a},y,{x,\[CurlyEpsilon],60}]; In[5]:= Plot[Evaluate[(x*y[x])^2+(((x-0.5 Pi)*y[x-0.5 Pi]))^2/.s[[1]]],{x,\[CurlyEpsilon]+Pi/2,60},PlotRange->All] (* see: http://dl.dropbox.com/u/3030567/Mathematica/schroed.png *) In[6]:= Evaluate[(x*y[x])^2+(((x-0.5 Pi)*y[x-0.5 Pi]))^2/.s[[1]]]/.x->{10,60} Out[6]= {0.0000979763,0.0000976444} does that help? Peter
From: Bob Hanlon on 3 Jul 2010 08:20 With PlotRange->All you are zoomed in and looking at numerical noise. Use PlotRange->{0,10^-4} to see that the result is essentially constant. Bob Hanlon ---- "ntina.86(a)libero.it" <ntina.86(a)libero.it> wrote: ============= Hello, I haven't been using Mathematica for some years and now I need to ask you for help with NDSolve. I'm using version 6.0. I want to solve the radial Schrodinger equation with a Yukawa potential, and I started trying to solve this: a = 500 b = 10; s = NDSolve[{y''[x] + (2/x) y'[x] + (1 + (2*a/x)*E^(-b*x)) (y[x]) == 0, y [0.000001] == 1, y'[0.000001] == -a}, y, {x, 0.000001, 60},PrecisionGoal -> \ [Infinity], MaxSteps -> Infinity] This gives me results that might be correct, but the problem comes when I try to plot: Plot[Evaluate[(x*y[x])^2 + (((x - 0.5 Pi)*y[x - 0.5 Pi]))^2 /.s[[1]]], {x, 10, 60}, PlotRange -> All] It should be a constant (i.e. the sum of the square of a cosine and of a sine, which is the solution of the Schrod.) at least for x>10, instead I get a badly oscillating function. Is it a problem of the method of integration or I'm doing something wrong? I tried to change it, using ExplicitRK or the SymplecticPartitionedRK, but the latter gives me these errors: --------------------------------------------------------------------------------------------------------- NDSolve::nlnum1: The function value {-21065.3 (1.23594*10^7+82.8327 \ TemporaryVariable[2][1.1462*10^-6])} is not a list of numbers with \ dimensions {1} when the arguments are {1.1462*10^-6,0.494377}. InterpolatingFunction::dmval: Input value {28.4292} lies outside the \ range of data in the interpolating function. Extrapolation will be \ used. >> InterpolatingFunction::dmval: Input value {30} lies outside the range \ of data in the interpolating function. Extrapolation will be used. >> NDSolve::nlnum1: The function value {-21193.5 (1.23496*10^7+82.3232 \ TemporaryVariable[2][1.14631*10^-6])} is not a list of numbers with \ dimensions {1} when the arguments are {1.14631*10^-6,0.493983}. InterpolatingFunction::dmval: Input value {28.4292} lies outside the \ range of data in the interpolating function. Extrapolation will be \ used. >> General::stop: Further output of InterpolatingFunction::dmval will be \ suppressed during this calculation. >> NDSolve::nlnum1: The function value {-21321.9 (1.23397*10^7+81.8195 \ TemporaryVariable[2][1.14643*10^-6])} is not a list of numbers with \ dimensions {1} when the arguments are {1.14643*10^-6,0.493589}. General::stop: Further output of NDSolve::nlnum1 will be suppressed \ during this calculation. >> ------------------------------------------------------------------------------------------------------- Thank you in advance. Regards, Valentina
From: Narasimham on 4 Jul 2010 03:10 On Jul 2, 4:27 pm, "ntina...(a)libero.it" <ntina...(a)libero.it> wrote: > Hello, > I haven't been using Mathematica for some years and now I need to ask you= for > help with NDSolve. I'm using version 6.0. --- > Thank you in advance. > Regards, > Valentina Boundary Conditions for y and y' need to be given at x = 0 exactly. If nearby points are taken for defining BC, different constants would result in. In the following, the constant looks wavy, but note y axis plot values are in a narrow range, are reasonably constant. s=NDSolve[{y''[x]+(2/x) y'[x]+(1+(2*a/x)*E^(-b*x)) (y[x])==0,y[0.1]==1,y'[0.1]==-a},y,{x,0.1,60}] Plot[Evaluate[(x*y[x])^2+(((x-0.5 Pi)*y[x-0.5 Pi]))^2/.s[[1]]],{x, 10,60},PlotRange->All] Narasimham
|
Pages: 1 Prev: non rectangular crop Next: How to edit a help reference page (e.g. NIntegrate.nb)? |