Re: Applying lists to FindRoot of a NIntegrate function
- To: mathgroup at smc.vnet.net
- Subject: [mg101381] Re: Applying lists to FindRoot of a NIntegrate function
- From: Alexei Boulbitch <Alexei.Boulbitch at iee.lu>
- Date: Sat, 4 Jul 2009 06:42:22 -0400 (EDT)
Dear members of the mathgroup email group,
I have serious difficulties in solving the following problem in Mathematica
(v.6) and would very appreciate it if anyone could give me a helping hand with
that.
To start with, I have a list of data points {{x[[i]],y[[i]]}}.
Also, I have the function W. W consists of two single functions f[r] and g[r,
x, a]:
f[r_] = r^2+r+2/(r^2+r)
g[r_, x_, a_] = Exp[a*Exp[-x*r]]
W := NIntegrate[f[r]*g[r, x, a], {r, 0.1, 100}]
What I have given for every pair of data is W[x[[i]]] == y[[i]] (see the list
of data points). What I do not have is the corresponding value of a[[i]].
So, I want to use FindRoot to give me the value of a[[i]] for each of my data
pairs {x[[i]],y[[i]]}.
Although I have found many examples quite similar to this problem, I have been
unable to adjust those to work in my case.
I would be very grateful if somebody could provide me with the proper syntax
for this.
Thanks very much in advance,
Chris
Dear Chris,
First, you have written your functions in a wrong way. Functions are fixed with ":=", rather than "=".
Second, you cannot use NIntegrate to integrate numerically a function with two unknown parameters.
Third, your expression W[x[[i]]] == y[[i]] namely, has sense, but not in your context.
Fourth, you may fix a nested list in the way you do {{x[[i]],y[[i]]}}, but typically it is better to do it differently.
All these things show that you did not do your homework. You should look more assiduously into your Mathematica
book or into tutorials before asking people on the forum.
Now to the point. I again have to repeat my favorite sentence: the problem here is not in Mathematica, but in mathematics.
Sorry. Before you understand what you need to do from the mathematical point of view, you should not start with Mathematica.
As much as I can understand you, what you have here is an equation, and you need to determine an unknown variable "a" from it.
The equation you get as the result of an integration, and a is incorporated into the integral in a complex way. Further, the
integral cannot be calculated exactly. (At least I made a fast trial with Mathematica and failed). To calculate the
integral one needs to substitute there a parameter x.
OK, the problem is that you cannot calculate the integral exactly, but on the other hand you cannot do it easily numerically
using NIntegrate, since it is nonintrinsic.
I can see here two approaches:
1. A more straightforward one is the following. First you should have a guess about the range in which a can lie.
Then you may make a table of different a values of a within this range with a fixed value of x. Then you get a set of values of
the integral. You fit it with some simple function (like a polynomial). Ans finally substitute this polynomial into equation
and solve it. This should be however, typically done for each x value and a interval separately, since you may have different
behavior of the integral, and also since your equations not necessarily have solutions within your expected intervals for a.
One example with a within the interval [0, 2], x=1 and y=338350 is below.
Before executing the code below remove a semicolumn after the Show operator. Your solution is a=1.47.
In[99]:= Clear[g, f];
f[r_] := r^2 + r + 2/(r^2 + r);
g[r_, x_, a_] := Exp[a*Exp[-x*r]];
tab = Table[{a, NIntegrate[f[r]*g[r, 1, a], {r, 0.1, 100}]}, {a, 0, 2,
0.1}]
(* this is the table of numerically calculated integrals *)
lp = ListPlot[tab];(* this is its plot*)
ft = Fit[tab, {1, a, a^2}, a] (*this is its fitting polynomial *)
pl = Plot[ft, {a, 0, 2}];(*this is the graph of the fit*)
Show[{lp, pl}];(* this shows the plot of the table of integrals and its fit \
together.Remove the semicolumn before the executing it *)
NSolve[338338 + 4.42 a + 2.53 a^2 ==
338350, a] (*this is the solution*)
Out[102]= {{0., 338338.}, {0.1, 338339.}, {0.2, 338339.}, {0.3,
338340.}, {0.4, 338341.}, {0.5, 338341.}, {0.6, 338342.}, {0.7,
338343.}, {0.8, 338344.}, {0.9, 338344.}, {1., 338345.}, {1.1,
338346.}, {1.2, 338347.}, {1.3, 338348.}, {1.4, 338349.}, {1.5,
338350.}, {1.6, 338352.}, {1.7, 338353.}, {1.8, 338354.}, {1.9,
338356.}, {2., 338358.}}
Out[103]= 338338.+ 4.42298 a + 2.53297 a^2
Out[104]= {{a -> -3.22003}, {a -> 1.47299}}
2. There is another way which I would chose, if it would be my research. Notice that the function g varies fast within a certain
interval 0<r<x. Outside of this interval it is mainly unity. I would try to analytically find an asymptotic approximation for this
integral basing on this fact. This analytical approximation I would substitute into
the equation and solve it then numerically. The method of the numeric solution should be chosen not earlier then the asymptotic
approximation for the integral is found.
Have fun, Alexei
--
Alexei Boulbitch, Dr., habil.
Senior Scientist
IEE S.A.
ZAE Weiergewan
11, rue Edmond Reuter
L-5326 Contern
Luxembourg
Phone: +352 2454 2566
Fax: +352 2454 3566
Website: www.iee.lu
This e-mail may contain trade secrets or privileged, undisclosed or otherwise confidential information. If you are not the intended recipient and have received this e-mail in error, you are hereby notified that any review, copying or distribution of it is strictly prohibited. Please inform us immediately and destroy the original transmittal from your system. Thank you for your co-operation.