First, is the curve fitting toolbox the best way to approach the problem? You COULD use it. Sadly, not really a good idea. That is because sums of exponentials are always difficult things to try to estimate, because the linear algebra very quickly encounters singularities. We can see why because one exponential, say exp(-x) has very much the same shape as another, for example exp(-2*x). So when we try to form sums of such terms, it becomes very difficult to extract the necessary coefficients using linear algebra. Next, you can only use fit on sets of points NOT as a way of approximating functions. So you could use hundreds or thousands of points on some interval, but that is all you could do. And they you are only finding a fit for a curve that passes through the listed set of points. Instead, I will want to approach such a problem using symbolic tools and function norms, etc.
A similar problem of near singularities happens in matrices like the Hilbert matrix, where we see such a matrix has very large condition numbers, even for relatively small matrices.
format rat
A = hilb(6)
cond(A)
plot(1:6,A,'-')
As you can see, even for a small matrix, the condition number is blowing up quite rapidly, even though the matrix is mathematically non-singular. But if you look at every row or column of the matrix, they ALL have very much the same shape. Each row looks almost like a scaled version of the other rows. And that is just a recipe for numerical problems when you venture into linear algebra. Sums of exponentials are similar.
Next, I never stated the domain of this approximation, and that is a CRUCIAL flaw in the problem statement. (The lack of a stated domain was intentional on my part.) Why? Suppose we decided to try to approximate the function 1/x in the entire first quadrant, thus from 0 to inf for x? The problem is, 1/x has a singularity at x==0, where it approaches infinity, but no such exponential will ever approach infinity at x==0. And no sum of such exponentials can do so either, at least with finite coefficients in front of each term.
In turn, that means we cannot do this approximation on the interval [0,inf]. Better if we try to do so on a more limited interval. Perhaps [1,inf] will suffice. It may have a chance of success. After all, the functions 1/x and exp(-x) do look sort of similar. At least at first glance, they do. So surely we can approximate one as a series of terms from the exponential family. No?
First, I'll see if I can even formulate and solve the problem for one term. This is always a good idea. Start small. Is there some fundamental problem in your approach?
syms a1 x real
approx = a1*exp(-x);
obj = int((1/x - approx).^2,x,[1,inf])
I formulated my approximation as a simple one term exponential, where only the unknown a1 is to be determined. The objective is the function norm of the difference between the two functions, so the integral of the square of the difference, on the half infinite interval. In fact, what I created in obj is known in mathematics as a functional, thus a function of an entire family of functions. I can now differentiate obj with respect to the parameter a1, setting that result to zero to determine the value of a1 which minimizes the functional.
First, does the integral even exist? As you can see, it is well defined, resolved by the function ei, which in MATLAB is the exponential integral, a special function. (For those of you who have a copy of texts like Abramowitz and Stegun.)
help ei
Now we can easily solve obj for the value of a1 which minimizes it.
a1sol = solve(diff(obj,a1) == 0)
vpa(a1sol)
Which makes sense, in that a1 should be approximately 1/exp(1). Plot the result, and see how we did.
fplot(subs(approx,a1,a1sol),[1,10])
hold on
fplot(1./x,[1,10])
hold off
And that does not look overtly terrible. After all, we only used one term in the sum. This was a VERY HIGHLY truncated series of exponentials. It did ok for small x, maybe below x==3.
I'll be more aggressive now though. Can we fit a higher order series, with perhaps 2, or 10 terms? Surely that will do better. More terms must be a good idea. (Ok, that is not always a good idea, if you read many of my answers here on this forum.) First, try 2 terms.
n = 2; % 2 terms in the approximation
a = sym('a',[1,n])
syms x real
approx = dot(a,exp(-(1:n)*x));
obj = int((1/x - approx).^2,x,[1,inf])
asol = solve(gradient(obj) == 0)
Again, we can see the solution uses ei extensively. A great little function when you need it. But that is the case for many of the special functions in mathematics.
fplot(subs(approx,asol),[1,10])
hold on
fplot(1./x,[1,10])
hold off
And while that looks a little better, it is clearly not that good. Even so, the general idea seems to work. Where 2 is ok, 10 MUST be better. (Will I ever learn?)
n = 10; % 10 terms in the approximation
a = sym('a',[1,n]);
syms x real
domain = [1,inf];
approx = dot(a,exp(-(1:n)*x));
obj = int((1/x - approx).^2,x,domain);
asol = solve(gradient(obj) == 0);
fplot(subs(approx,asol),[1,10])
hold on
fplot(1./x,[1,10])
hold off
The result has the essential character of what such a curve fit should look like. The two functions will cross back and forth. If I plotted the difference, we would see an oscillating behavior back and forth around zero, as least for small x. Beyond 10, the oscillations will stop, as the difference very slowly dies off towards zero.
fplot(1./x - subs(approx,asol),[1,10])
yline(0)
But do you see the problem? Terms of the general form exp(-k*x) just go to zero so much faster than 1./x. While we can get arbitrarily good approximations for reasonably small limits on x, at some point, all of those terms just die off. For significantly larger values of x, we might need hundreds or even thousands of terms in that approximation.
If we truly needed a good approximation with a finite number of terms in the series, we need to limit the domain. If not, this is not any different from trying to use a Taylor series approximation for functions like sin(x) or exp(x). While they are globlly convergent in theory as long as you take enough terms, we will never be able to take as many terms as we really need in practice for sufficiently large x.
Have I now completely killed this problem? Perhaps. The next extra step I might consider is to perform the approximation on a finite domain, perhaps [1,5].
n = 10; % 10 terms in the approximation
a = sym('a',[1,n]);
syms x real
domain = [1,5];
approx = dot(a,exp(-(1:n)*x));
obj = int((1/x - approx).^2,x,domain); % function norm on a finite interval
asol = solve(gradient(obj) == 0);
fplot(subs(approx,asol),domain)
hold on
fplot(1./x,domain)
hold off
Wow! That was actually surprisingly good. When I first did this example, I restricted the domain to [1,3], and with 10 terms the fit was so good I could not see the difference between the two functions in the plot. So I expanded the domain to go as high as 5. Even on the larger domain, the two can be seen to be reasonable approximations to each other.
fplot(1./x - subs(approx,asol),domain)
yline(0)
And here we see the characteristic osculatory error plot for such an approximation, bouncing back and forth across zero.
Is this finally time to have killed off the question in my eyes? SIGH. Probably not. With some effort, I might even be able to derive a general expression for the nth such term in the series, as an infinite series. I'd not be surprised if such a solution does exist. In fact, a quick search finds this:
But for now, even though that might be fun, I've gone far enough.
