# Non-integrable multiple integrals for Mathematica

• Mathematica
“Non-integrable” multiple integrals for Mathematica

Dear all,

I have been trying to crack one problem in Mathematica, but I keep getting a wrong answer probably because I have something either fundamentally wrong analytically or code wise. OK, here is the problem.

Suppose we have to evaluate the following:

Integral from 0 to t, dm [of Integral from m to t, dd [of Integral from d to b of (y-d), dy],

and suppose we know that b = 3. Clearly, it is straightforward to accomplish this task (please see the first 3 steps in the code attached (Example2.nb files) and the “final” answer). (I deliberately proceeded in simple steps to keep a better track of the individual answers).

Suppose also I know that t is a time variable which varies from, say, 0 to 2. Then I can easily take various values of t from that interval, plug them into the “final” answer and plot if I need to. (Step size 0.01 would be precise enough for my purpose).

However, suppose we want to plot the “final” for "all" possible values of t, without being able to analytically integrate the first inner integral, i.e. “Integral from d to b of (y-d), dy”. Let us just pretend that the reason is that this thing is “non-integrable”.

I have performed a sequence of simple steps to get to the “finalnew”, but you can see that this is a wrong array, because if I set t=0.01, for example, in “final” I would get 0.000224001, which is not what the 2nd element in “finalnew” array is.

I would really appreciate if you could help. This is for my research, where the primary problem is that in place of (y-d) I have something pretty-complicated, which is “non-integrable”. I can numerically integrate it given that d was a dummy for t which goes from 0 to 2 as you recall. But then the question is, given 2 remaining integrals, how to build a final array of numbers to plot. All I need is the final array of right numbers. I am sure there is a correct and perhaps an easier way to accomplish this task, but the question is what it is.

Thanks a lot for your time.

Anna.
-----------------------------------------------

#### Attachments

• 11.6 KB Views: 368

Related MATLAB, Maple, Mathematica, LaTeX News on Phys.org

Would it be possible for you to make the second half of your notebook significantly less complicated?

I have made a serious effort to go through your notebook three times, each time starting from scratch trying to decompile it back into something much simpler that I can understand and thus figure out where the error is. And I have failed.

You do not have any outright syntax errors, which is better than many questions posed here, but you have integrals of lists of constants that you transpose into matricies only then to turn back into vectors by applying parts of the row to the rest of the row.

Your example in the first half, where you can use Integrate, simplifies down to Integrate[y - d, {m, 0, t}, {d, m, t}, {y, d, b}]

For the second half it is clear you have spent a lot of work trying to find a way to force the available Mathematica tools into a form that will try to calculate what you want. But is there any more conventional Mathematica construct you could use instead of repeatedly using Apply on the Transpose on the Integrate. I'm thinking this must represent something simpler, perhaps possible to create with something as simple as a Table[NIntegrate[something],{range},{range}].

If not then perhaps a simple English description of what you are trying to accomplish with each step might make it feasible to find a simpler and more easily understood way of doing this that would have as a side effect also being correct.

Or as a last resort, could you just manually calculate and show exactly what the values in the first half dozen or dozen entries of finalnew should be and precisely how you calculate the corresponding different values in final?

At one time I had almost reached the conclusion that part of the problem might be that your first NIntegrate always uses b as the upper limit of integration, rather than integrating only over an interval of width step, but after another attempt from scratch to decompile what you are really trying to accomplish I am not even convinced that is correct.

Last edited:

Hi Bill,

I have, BTW, just figured out what my error was. I failed to take into the account the fact that the first step should represent not just a single number, but the multiple areas under the curve. My approach was mathematically wrong.
Unfortunately, that confused me even more, frankly, i.e., now I am unsure how to crack the problem.
I have attached a pdf file, which explains the problem, and an nb file, where I have just assumed that F(y,d) is y-d for the sake of illustration. All I need is get to the answer in the form of a “finalarray” as shown, but without using Integrate (assuming that y-d cannot be Integrated directly).
I have tried using interpolating function, etc. but my poor programming skills result in being stuck.

Anna.
P.S. This is NOT a hw assignment, but this just would greatly simplify my research as I deal with this type of problem a lot.
P.P.S. I have deliberately didn't insert into an nb file an original function F(y,d) since I know for sure it is "non-integrable", and it would occupy like half a page. If I know the procedure for a simple y-d, then I can apply it for other functions.

#### Attachments

• 32.2 KB Views: 160
• 13.9 KB Views: 398

You ask for the final answer in the form of an array. Is that array containing the integral for various values of the unknown t? The unknown m? (I'm guessing not) I'm trying to clearly identify exactly what the index of that array represents, understand exactly the form of the "prism" that you are integrating over and clearly see how your function F relates to the individual dVolume bits making up the prism.

Not showing pages of meaningless detail for your F(y,d) is perfectly fine, better than confusing the issues by showing that. If, for the moment, you just do not define f[y,t] you could insert that into your Mathematica code and not try to actually evaluate the result, but just show what you are hoping to calculate.

NIntegrate, like the other Nfunctions... NSolve and NDSolve, is strict in its' demands that everything evaluate to numeric constants. You cannot expect NIntegrate to evaluate any expression that contains free variables.

I assume, and please confirm this, that f[y,d] is not independent of y and d, otherwise why would you have parameters y and d, but is independent of m and t. Thus, in light of the previous paragraph, I believe you cannot use NIntegrate on f[y,d] for a single value of y and d and then use that result across all values of y and d. But if my assumption is correct then you can NIntegrate your innermost f[y,d], that would result in a scalar constant and the two outer integrals appear to be over a triangular area, so you multiply that scalar times the area of the triangle and you would be done.

But your second level integration involves dd and I can't quite see yet how that is going to interact with your innermost integration.

Can you explain some of the things I still do not understand?

Last edited:

Hi Bill,

1. Yes, you are right, I need a final array of numbers for various values t. From the pdf file if we assume that F(y,d) is y-d, then (1) results in 9t^2/4 –t^3 + t^4/8. I do not need to “see” what the final explicit answer is in terms of a mathematical function of t, but I do need to know the “approximate” values of the function for the values of t from 0 to 2 in step size 0.01.
2. Yes, f(y,d) depends only on y and d, and not on m, t or anyhting else.
3. By assumption, b is given and d takes trhe values from 0 to t(=2) and F(y,d) is such that it is possible to aplly NIntegrate.

I think my problem is the coding issue. That is, given the above assumptions, how to code the problem which would construct that final array of 201 elements, without applying Integrate command.
Thanks.
Anna.

Just out of curiosity, exactly what does the following give you when you evaluate it?

f[y_, d_] := yourlargefunctiondefinition;
Table[NIntegrate[f[y, d], {m, 0., t}, {d, m, t}, {y, d, 3.}], {t, 0., 2., .25}]

Last edited:

Dear Bill,
Thank you very much for your reply. Your method is indeed very efficient and clever, and I have got the results, all being real numbers, with no hassle in the computing NIntegrate process...I appreciate the time you have invested into this question.
I have tried to extend your method to a “slightly” different problem, which I also often face, but perhaps because I making a logical/procedural error, I get stuck.
Can you please see the attached pdf file? Essentially, the problem “looks” similar to the previous one, but this time the inner integral is an exponent of a double integral. All other assumptions remain the same. I am attaching the nb file too.
I tried various combinations, and even Log substitution for an inner integral, but nothing helps...Perhaps I need to develop a better understading of how Mathematica handles multiple integrals, but so far I didn't come accross to a problem very similar to this one ...
Thanks a lot,
Anna.

#### Attachments

• 10.7 KB Views: 378
• 33.4 KB Views: 149

I'm glad we were able to find such a simple solution to your initial problem. Just to be cautious, you should carefully check the numerical results of my NIntegrate expression to make sure that it is giving you right answers.

Last edited:
Dale
Mentor

NIntegrate is definitely the way to go, and it is much better to do it as a single command than to split it up into three integrals. The numerical routines can control the overall accuracy and precision if you evaluate it as a single NIntegrate. Be sure to read the documentation on NIntegrate so that you understand what it is doing.

Hi Bill,
Just to make sure that I understood you correctly, you are asking me to verify if (1) (please see the attachment) is true or not? I will be working on it to see if I can verify it.
Thanks a lot.
Anna.

#### Attachments

• 20.8 KB Views: 138

I have made an error in what I was trying to accomplish.

Last edited:

Now as an alternate formulation multiply each dy*dd*dm by E^m, integrate over dy and dd to get E^m times the area and finally integrate over dm to get your weighted volume.

Hi Bill,
Sorry for a silly question, - I have been trying various techniques to get this problem into a simpler looking triple integral according to your suggestions, - so I might be a bit slow ... When you say "multiply each dy*dd*dm by E^m", do you mean that you "insert" E^m into each of the 3 integrands? And then you use the integration by parts for every innermost integral?
Thanks a lot.
Anna.

Bill,
Correct me if I am wrong, but it appears to me that when we try to make E(of double integral) become of the form of a double integral (of some function), then essentially weare solving an integral equation. If so, I fear that this might not be an easy task...But I might be wrong ...
Thanks a lot.
Anna.

Ah, no. I am sorry but I believe I have led you down a path that does not work. I think I see now that the idea is incorrect and cannot be repaired.

I am very sorry for my error and taking up your time with that.

So I think you are back to your Integrate[E^Integrate[]] problem. If you can try to show exactly how you are having a difficulty evaluating that then perhaps we can find a solution for you.

Last edited:

Ah, no. I am sorry but I believe I have led you down a path that does not work. I think I see now that the idea is incorrect and cannot be repaired.

I am very sorry for my error and taking up your time with that.

So I think you are back to your Integrate[E^Integrate[]] problem. If you can try to show exactly how you are having a difficulty evaluating that then perhaps we can find a solution for you.
Dear Bill,
Please do not apologize as I really appreciate your time and very useful hints.
I will try to formulate a solution, and see if it can get us somewhere...I will be in touch soon. My only concern now is as follows: since the problem is back to the original, i.e. cannot be represented as a triple integral, do you think I must use some kind of nested looping to get to the answer, or you beleive a proper representation of NIntegrate[etc.etc...] can still work?
Many thanks.
Anna.

Dale
Mentor

This should work provided everything is well behaved:

NIntegrate[Exp[NIntegrate[f[y, d], {d, m, t}, {y, d, b}]], {m, 0, t}] /. {f -> Function[{y, d}, y], t -> 5, b -> 2}

Just set your particular values for f, t, and b. It will throw some warnings about non-numerical values which you can ignore.

Last edited:

This should work provided everything is well behaved:

NIntegrate[Exp[NIntegrate[f[y, d], {d, m, t}, {y, d, b}]], {m, 0, t}] /. {f -> Function[{y, d}, y], t -> 5, b -> 2}

Just set your particular values for f, t, and b. It will throw some warnings about non-numerical values which you can ignore.
Dear DaleSpam,
Thanks you very much for your response.
In my problem (please see the attachment) t is time variable, which is assumed to go from 0 to 2 (step size in 0.01) is sufficient for the accuracy.
That is the final answer is an array of 201 elements. I wonder if your method can be used where t goes from 0 to 2 in step 0.01. Maybe combining it with a Table command?
I am playing with it now to see how far I can get...
Thanks a lot.
Anna.

#### Attachments

• 28.7 KB Views: 146

P.S. Hi DaleSpam,
I have tried this Table[NIntegrate[
Exp[NIntegrate[f[y, d], {d, m, t}, {y, d, b}]], {m, 0, t}] /. {f ->
Function[{y, d}, y], b -> 3}, {t, 0, 2, 1/4}], but I fear I might be wrong...

Dale
Mentor

That ran on my computer.

One other thing to think about. If your function is fairly smooth then you may have a computational advantage by making an interpolating function for the argument to the Exp function. However, I am not certain of that.

Last edited:

This should work provided everything is well behaved:

NIntegrate[Exp[NIntegrate[f[y, d], {d, m, t}, {y, d, b}]], {m, 0, t}] /. {f -> Function[{y, d}, y], t -> 5, b -> 2}

Just set your particular values for f, t, and b. It will throw some warnings about non-numerical values which you can ignore.
Hi DaleSpam,
Sorry for a silly question, but what does {f -> Function[{y, d}, y] exactly accomplish in the ab ove code? It appears that without the task can be accomplished too, but maybe I misunderstand something...Can you please also advice where I can read more about that syntax?
Many thanks.
Anna.

Dear DaleSpam and Bill,
Thanks a lot for all your responses as they have been extremely helpful and made me learn a lot. I appreciate both of yours time.

To make sure that I have a proper understanding, I started twisting the problem a little bit (and also because exp problem is triky to symbolically integrate), and it seems to me that the method you have offered works.

If you see the pdf file and nb file, and it looks like it works fine and it given me an error (please see “diff” value less than 1/1000000).

One thing I am worried about is the property of F being well-behaved, as DaleSpam has indicated that the method would work in such a case. I wonder if the method would not for more complex functions...

1) why so I get those “\[NoBreak]y\[NoBreak] = \[NoBreak]d\[NoBreak] is not a valid limit of integration. ” errors and if they pose any “danger” really? I also wonder why despite those errors the right answer is reached anyway.

2) Further to my very last post, as I have deleted f-> function ... syntax, I wonder if this is OK to do. Maybe yes for this particular one, but not the one where had this exponent of double integral?

Thanks a lot.
Anna.

#### Attachments

• 28.1 KB Views: 130
• 12 KB Views: 410
Dale
Mentor

Sorry for a silly question, but what does {f -> Function[{y, d}, y] exactly accomplish in the ab ove code?
That is just where the function f(y,d) is defined. For testing purposes I just had f(y,d)=y.

Dale
Mentor

1) why so I get those “\[NoBreak]y\[NoBreak] = \[NoBreak]d\[NoBreak] is not a valid limit of integration. ” errors and if they pose any “danger” really? I also wonder why despite those errors the right answer is reached anyway.
Those are due to the order of evaluation. First it tries to evaluate the NIntegrate, it finds that there are still some symbolic limits, and throws the error. Then it makes the substitution, re-evaluates NIntegrate, finds that there are no symbolic limits and succeeds.

Hi DaleSpam,

It seems the method you have suggested is really flexible and efficient and can be applied to a number of variations succesfully!
Unfortunately, as those variations get messy enough, it becomes impossible for me to re-check whether the answer is exactly the same as “Integrate” commands would imply (especially with Exp).
I have just wrote a small code (please see the nb attachment) for the problem in the pdf file. However, I fear that maybe I have made a mistake due to a potential misunderstanding of your idea. Can you please have a quick look to see if what I wrote is indeed a correct representation of the problems solution, and does not go against your idea? The code obviously runs well, but I fear that I have made an error.
Anna.

#### Attachments

• 34.8 KB Views: 128
• 10.7 KB Views: 442
Dale
Mentor

It looks pretty good. Here are a couple of small suggestions.

Because you are repeating that inner-most integral so many times it is actually quite a bit more efficient to evaluate it once, construct an interpolating function and then just evaluate that interpolating function in the other integrals. It calculates about 26 times faster with little loss in accuracy.

#### Attachments

• 8.2 KB Views: 417