Group Abstract Group Abstract

Message Boards Message Boards

0
|
6K Views
|
21 Replies
|
6 Total Likes
View groups...
Share
Share this post:

Accurately solving an ode/bvp?

Posted 9 years ago

Hello,

I have been struggling (actually for several years) to accurately solve a BVP similar to the following BVP, using MATHEMATICA:

y"(x) - (s + x)*y(x) == 0,  y'(0) == -1, y(infinity) == 0

, where s is a positive parameter from the range between 0 and 10^2 (let's say) or higher. The above example possesses an analytical solution in terms of Airy functions, but my problem (slightly different) can only be solved numerically. I need a very high accuracy (relative errors 10^(-32) or smaller), but MATHEMATICA seems to fail to achieve even much worse accuracy, despite the propaganda that it can in principle provide an arbitrary accuracy of any calculations. The problem formulated as a BVP (also with infinity replaced by a finite distance) apparently cannot be solved with errors smaller than about 10^(-n) where n = 2-5 or so. Somewhat better effects can be achieved by converting the BVP to an initial value problem and using the shooting method (programmed by hand), but the accuracy is still insufficient. I have not been able to obtain correct results by using the shooting method built into the BVP solver of recent versions of MATHEMATICA. So, I definitely need an expert advice how to proceed in order obtain the relative errors 10^(-32).

Leslaw

POSTED BY: Leslaw Bieniasz
21 Replies

BTW, the first two posts where I just repeat your last post were not posted by me.. I think. I just posted the last one and was close to taking a video of me submitting my post and getting reposts of yours out of it.

Perhaps the moderation team could delete the last two re-posts of Sanders post. They were not intended, and I have no idea how they got there.

Cheers,

Marco

POSTED BY: Marco Thiel

I'll reply to Sander's post before I can see it:

Yes, I have read that, too. Would be good to know how this works exactly. I've got several computers on my desk. The all update slowly, but at different times.

Also, I am not complaining, but we are currently using this thread nearly like a chat. I guess that's why the delay becomes more of an issue. In understand that the Community is not quite a chat-room, which might be another idea for collaborative work. ;-)

Cheers,

Marco

POSTED BY: Marco Thiel
POSTED BY: Leslaw Bieniasz
POSTED BY: Marco Thiel

Have you read:

https://reference.wolfram.com/language/tutorial/ControllingThePrecisionOfResults.html

? Notice the phrase:

"But within this constraint you can tell the Wolfram Language how 
much precision and accuracy you want it to try to get."

In order to attain a certain precision/accuracy you need a certain workingprecision and a good algorithm to 'home in' to the correct solution. So two things are possible:

Chose an algorithm explicitly which improves the way it 'homes in' to a good solution. A crappy integration scheme limits the possible accuracy/precision one can get when given a certain working precision. By default it accuracy/precision goal to be half the digits of the workingprecision. If accuracy/precision goals becomes closer to workingprecision you need smarter and smarter integration schemes in order to get the correct result. It is sometimes much faster 'just' increasing the workingprecision, and set no goals (defaults) rather than trying different algorithms/integration schemes and stepsizes and so on...

POSTED BY: Sander Huisman

I experience huge delays on the update of the community. I receive emails of your posts and they only show up in the browser long after that.

Cheers,

M.

POSTED BY: Marco Thiel
POSTED BY: Frank Kampas

To paraphrase the creators: the site has multiple servers/databases which are synchronised slowly. I see yours instantly, but that just means I'm logged on to a 'good one' I guess. Sometimes it is also slow for me. Not this time though... Speed improvements would be welcomed indeed @Vitaliy Kaurov. Especially the groups and the main overview load very slow.

POSTED BY: Sander Huisman

Hi everyone,

It appears that there is no issue with the analytical solution here, right? (I might have overhead something, because I could not read all posts, as I am on a mobile phone.)

sol = DSolveValue[{y''[x] - (s + x)*y[x] == 0, y'[0] == -1, y[M] == 0}, y[x], x]
(*(-AiryAi[s + x] AiryBi[M + s] + AiryAi[M + s] AiryBi[s + x])/(AiryAiPrime[s] AiryBi[M + s] - AiryAi[M + s] AiryBiPrime[s])*)

Then we can calculate

limit=Limit[sol, M -> Infinity]

enter image description here

It seems to be ok-ish, from what I can see now:

D[limit, x, x] - (s + x)*limit == 0 

evaluates to True.

D[limit, x] /. x -> 0

gives -1; and

Limit[limit, x -> Infinity]

gives 0. From that we can generate very high precision results, of course.

But as the OP says, in the case they are interested in there is no analytical solution. So in order to suggest anything it would be useful to get the actual problem.

By the way, the example:

sol=NDSolve[{u''[x]-(s+x)*u[x]==0,u'[0]+1==0,u[L]==0},u[x],{x,0,L},WorkingPrecision->50,PrecisionGoal->32];

fails to run on my MMA11.0.0 on OSX 10.11.6, but it does run on MMA 10.4.1. I also get an error message that the maximum number of steps is reached before obtaining the requested precision, but that can be changed by adding something like this:

post = NDSolveValue[{u''[x] - (s + x)*u[x] == 0, u'[0] + 1 == 0, u[L] == 0}, u[x], {x, 0, L}, WorkingPrecision -> 50, PrecisionGoal -> 32, MaxSteps -> 10^6]

which does finish after quite some while. No joy on MMA11. There are several reasons why I think that your way of checking the precision is not quite ideal, but this is what I get using your approach:

rerror = u0/u0acc - 1

gives $7.1796\cdot 10^{-27}$. Not quite what you want though.

Finally, if you use

solpost = 
 NDSolveValue[{u''[x] - (s + x)*u[x] == 0, u'[0] + 1 == 0, u[L] == 0},u[x], {x, 0, L}, WorkingPrecision -> 56, PrecisionGoal -> 39, 
  MaxSteps -> 10^6, InterpolationOrder -> All, Method -> "StiffnessSwitching" ]

your test gives basically what you might expect $3.335\cdot10^{-32}$. You can play a bit with the parameters to improve on this.

And this here

solpost = 
 NDSolveValue[{u''[x] - (s + x)*u[x] == 0, u'[0] + 1 == 0, u[L] == 0},u[x], {x, 0, L}, WorkingPrecision -> 60, PrecisionGoal -> 45, 
  MaxSteps -> 10^7, InterpolationOrder -> All, Method -> "StiffnessSwitching" ]

gives $1.1*10^{-34}$ for your "rerror".

Cheers,

Marco

POSTED BY: Marco Thiel

"Not quite what you want though."

Also note the AccuracyGoal documentation page:

"Even though you may specify AccuracyGoal->n, the results you get may sometimes have much less than n-digit accuracy."

AccuracyGoal is a goal, but it might not reach it (everywhere).

POSTED BY: Sander Huisman

Hi Sander,

yes, I had that clear. The point is that the first approach might be to crank up the Accuracy/Precision and hope that it'll work. But it did not. Needs e.g. explicit method for integration etc.

Cheers,

Marco

POSTED BY: Marco Thiel

Hello,

OK, you are right, increasing WorkingPrecision helps, so that I wrote a new version of my test code:

(* Analytical solution, assuming finite interval L, just to avoid any potential controversy regarding the right choice of L to mimic infinity *)

   uanal[x_,s_,b_]=-(AiryAi[s+x]*AiryBi[s+b]-AiryAi[s+b]*AiryBi[s+x])/(AiryAiPrime[s]*AiryBi[s+b]-AiryAi[s+b]*AiryBiPrime[s]);

   s=2;

   L=10;

   ndigits=32;

   sol=NDSolve[{u''[x]-(s+x)*u[x]==0,u'[0]+1==0,u[L]==0},u[x],{x,0,L},WorkingPrecision->50,PrecisionGoal->30];

   uappr[x_]=First[u[x]/.sol];

  u0=N[uappr[0],ndigits]

        0.65782402587273264492979577402005

  u0acc=N[uanal[0,s,L],ndigits]

       0.65782402587268160074874669507081

  rerror=u0/u0acc-1  

       7.75954952106876733e-14

As can be seen, a relative error about 8e-14 is now achieved. However, please note that I requested AccuracyGoal->30. I don't see a way to obtaining such accuracy, but I need at least 32 accurate digits. Setting AccuracyGoal->32 yields the following solution failure:

sol=NDSolve[{u''[x]-(s+x)*u[x]==0,u'[0]+1==0,u[L]==0},u[x],{x,0,L},WorkingPrecision->50,PrecisionGoal->32];

    NDSolve:mxst: Maximum number of 10000 steps reached at the point x== ...etc.
    NDSolve:bvaux: Unable to solve auxiliary system for boundary conditions at x == 0

So, Please tell me how to achieve my goals.

Leslaw

POSTED BY: Leslaw Bieniasz

If I run your code as-is:

uanal[x_, s_, b_] = -(AiryAi[s + x]*AiryBi[s + b] -  AiryAi[s + b]*AiryBi[s + x])/(AiryAiPrime[s]*AiryBi[s + b] - AiryAi[s + b]*AiryBiPrime[s]);

s = 2;
L = 10;
ndigits = 32;
sol = NDSolve[{u''[x] - (s + x)*u[x] == 0, u'[0] + 1 == 0, u[L] == 0}, u[x], {x, 0, L}, WorkingPrecision -> 50, PrecisionGoal -> 30];
uappr[x_] = First[u[x] /. sol];
u0 = N[uappr[0], ndigits]
u0acc = N[uanal[0, s, L], ndigits]
rerror = u0/u0acc - 1

I get:

-2.06657*10^-26

Using version 10.4.1 or version 11 same result.

If you want more steps you can set MaxSteps option. Look at the documentation! This is all quite trivial if you just spend a few minutes reading and observe the examples.

POSTED BY: Sander Huisman

Even the error messages have documentation:

ref/message/NDSolve/mxst

which explicitly tell you that you can increase MaxSteps.

Also note the AccuracyGoal documentation page:

"Even though you may specify AccuracyGoal->n, the results you get may sometimes have much less than n-digit accuracy."

POSTED BY: Sander Huisman

Hi Sander,

I find it a bit worrying that that code works on your machine but not on mine; doesn't work on MMA11.

Cheers,

M.

POSTED BY: Marco Thiel

What about 10.4(.1)? I get the same result in both.

POSTED BY: Sander Huisman

10.4.1 works fine.

Performed clean start on MMA11; still the same. Will investigate further.

M.

POSTED BY: Marco Thiel
POSTED BY: Leslaw Bieniasz

I don't know u0acc so I don't know what is the correct solution. But:

s = 10; L = 10; sol = 
 NDSolve[{u''[x] - (s + x)*u[x] == 0, u'[0] + 1 == 0, u[L] == 0}, u, {x, 0, L}, WorkingPrecision -> 50];
sol = u /. First[sol];
sol[0]

works just fine, without any error:

0.31379386687943471460767057577019210318851451370584
POSTED BY: Sander Huisman

Did you change the independent variable so that the range of the problem is from 0 to 1 rather than 0 to infinity?

POSTED BY: Frank Kampas

Without an example and the code you tried no one will be able to help you effectively, post your code, and your analysis of the n = 2--5.

POSTED BY: Sander Huisman
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard