Solving a system of integral equations numerically

I'm trying to find the numerical solution to a system, the first problem I have is that I'm given a system of two integrals with two unknowns where the unknowns are in the integrands and the limits of both equations. How do they do this? I've attached the problem below. I've tried just using the built-in solve function but it doesn't work for this integral as I'm guessing its too complicated. I'm not too great on MATLAB, so any help would be appreciated! Thanks

Answer 1

[Star Strider]

There are several roots, and it doesn’t always converge on the published solutions.

The Code —

a2 = @(b) integral(@(y) 1./sqrt(((y.^2 + b(2))./b(1)).^2 - 1), 0, sqrt(b(1)-b(2))) - 1;
b2 = @(b) integral(@(y) ((y.^2 + b(2))./b(1))./sqrt(((y.^2 + b(2))./b(1)).^2 - 1), 0, sqrt(b(1)-b(2))) - pi/2;

B = fsolve(@(b) [a2(b), b2(b)], rand(2,1))
ym = sqrt(B(1) - B(2))

The Solutions —

B =
 -1.034659029995483 + 0.000000635413460i
 -2.290004434974137 + 0.000000465384426i

ym =
  1.120421976301188 + 0.000000075877231i

The imaginary parts are vanishingly small, so you can likely just ignore them.

[Roger Stafford]

@Star: Producing complex values is presumably caused by using the wrong initial estimates, x0, leading to complex values in the problem's integrands or integration limits. I warned about this in my answer. You should not trust any answer in which complex values are produced.

[Star Strider]

The imaginary parts are on the order of 1E-7 of the real parts, so approximately sqrt(eps), and essentially 0. This is obviously due to floating-point approximation error. With extended precision, more constrained tolerances, and more iterations, the solution will likely produce an entirely real result.

[Lewis Hancox]

I'm trying to get your code to run Star, but it says "Failure in initial objective function evaluation. FSOLVE cannot continue." I've literally copied both of the boxes out. Do they need to be separate or? Thanks again

[Torsten]

Try

x0 = rand(2,1);
x0 = [max(x0);min(x0)];
B = fsolve(@(b) [a2(b), b2(b)], x0);

Best wishes

Torsten.

[Lewis Hancox]

Still the same error, this is what I'm putting in now.

if true
  % code
end
a2 = @(b) integral(@(y) 1./sqrt(((y.^2 + b(2))./b(1)).^2 - 1), 0, sqrt(b(1)-b(2))) - 1;
b2 = @(b) integral(@(y) ((y.^2 + b(2))./b(1))./sqrt(((y.^2 + b(2))./b(1)).^2 - 1), 0, sqrt(b(1)-b(2))) - pi/2;
x0 = rand(2,1);
x0 = [max(x0);min(x0)];
B = fsolve(@(b) [a2(b), b2(b)], x0);
ym = sqrt(B(1) - B(2))

[Star Strider]

@Torsten — Thank you!

Another option is:

B = fsolve(@(b) [a2(b), b2(b)], rand(2,1)*10)

It is usually necessary to run the code a few times, with random initial parameter estimates, in order for it to converge successfully. All nonlinear parameter estimation problems are sensitive to the initial parameter estimates. Derivative-free solvers (such as genetic algorithms, the Global Optimization Toolbox ga function) avoid these problems, at the expense of slower convergence.

[Lewis Hancox]

I've pressed run so many times now, and the same error keeps coming now. How did you get the solution so easily?

[Torsten]

Maybe

if true
  % code
end
a2 = @(b) integral(@(y) 1./sqrt((y.^2 + b(2))./(b(1)^4 + b(2)).^2 - 1), 0, b(1)^2) - 1;
b2 = @(b) integral(@(y) ((y.^2 + b(2))./(b(1)^4 + b(2)))./sqrt((y.^2 + b(2))./(b(1)^4 + b(2)).^2 - 1), 0, b(1)^2) - pi/2;
B = fsolve(@(b) [a2(b), b2(b)], rand(2,1));
ym = b(1)^2

works better.

Best wishes

Torsten.

[Star Strider]

I am running R2017b, so there could be version differences.

If you use the randi function (introduced in R2008b) to choose two negative integers randomly, it converges, and also eliminates the complex number problem:

B = fsolve(@(b) [a2(b), b2(b)], -randi(9, 2, 1))

[Lewis Hancox]

Finally got it to work! Thanks guys, I had two open with similar codes so they were fighting each other.

[Lewis Hancox]

Quick follow up now, now I have the integral

if true
  % code
end
x = 1./sqrt(((Y.^2 - 2.29)./-1.035).^2 - 1), 0, 1.12) - 1;
I would I go about plotting this, with the values we found (-2.29,-1.035,1.12) Thanks

[Star Strider]

My (our) pleasure!

I don’t see an integral call, so I don’t know what you’re integrating. I don’t understand what you want to do with this:

x = 1./sqrt(((Y.^2 - 2.29)./-1.035).^2 - 1), 0, 1.12) - 1;

Are you:

(1) substituting the solved values for ‘B’ and ‘ym’ for ‘Y’,

(2) solving the integral of something for ‘Y’, using the solved values for ‘B’ and ‘ym’ as ‘x’,

(3) or something else, such as supplying a vector of values for ‘Y’?

Plotting it is straightforward once you explain what you want to plot.

Note that the integral of one point is a single value so it wouldn’t be a very informative plot.

I would also use the best precision available to you, not simply three or four significant figures. I include them here for reference:

B = [-1.034658762541066
     -2.290003426524971];
ym = 1.120421645624496;

[Lewis Hancox]

Oh sorry yes, its

if true
  % code
end
x=int(1./sqrt(((Y.^2 - 2.29)./-1.035).^2 - 1), 0, ym) - 1)DY

[Star Strider]

Here you go:

ym = 1.120421645624496;

x = integral(@(Y) 1./sqrt(((Y.^2 - 2.290003426524971)./-1.034658762541066).^2 - 1), 0, ym) - 1

x =
     3.206812593248287e-11

[Lewis Hancox]

Oh, I think the top limit is y and that varies.. So then I could plot it. Like introduce a linspace y.

[Star Strider]

Try this:

ym = 1.120421645624496;
ymv = linspace(0, ym);
for k1 = 1:numel(ymv)
    x(k1) = integral(@(Y) 1./sqrt(((Y.^2 - 2.290003426524971)./-1.034658762541066).^2 - 1), 0, ymv(k1)) - 1;
end

figure(1)
plot(ymv, x)
grid

[Lewis Hancox]

It looks good, I've switched the x and ymv in the plot(ymv,x) and I get one half of what I'm looking for, but if I extend ymv past ym in the linespace it gets messy. When it should just be a mirror image. I've attached the plot I got and the plot I should be getting

[Star Strider]

I have no idea how the authors get that result.

The only way I can get it is to fudge it:

ym = 1.120421645624496;
ymv = linspace(0, ym);
for k1 = 1:numel(ymv)
    x(k1) = integral(@(Y) 1./sqrt(((Y.^2 - 2.290003426524971)./-1.034658762541066).^2 - 1), 0, ymv(k1));
end

figure(1)
plot(x, ymv,  (2*max(x)-x), ymv)

[Lewis Hancox]

Exactly what I was thinking, I just wanted to add another curve on it too. It seems impossible, the square root goes negative and then its imaginary for a bit, and then increasing again. Thanks for that anyway!

[Star Strider]

As always, my pleasure!

Answer 2

[Roger Stafford]

This is a problem you can solve using Matlab's 'fsolve' function. The fact that your objective function requires two numerical integrations at each evaluation does not change the nature of the problem. It only means that execution time will be longer than with simpler objective functions. Each iteration requires the numerical evaluation of two integrals. You can presumably use Matlab's 'integral' function in evaluating the objective function.

You may have to be careful as to the selection of your initial estimate, x0, since it is easy to run into complex values with your particular integrals, both in the integrand and the integral limits.

Answer 3

[Roger Stafford]

@Lewis: I noticed that the integrands in both of your integrals become infinite at the upper limit of integration, and this can result in some inaccurate numerical results from the integration process. However, with the appropriate change of variables these integrals can be converted to complete elliptic integrals of the first and second kinds for which Matlab has functions that you can call on directly and thereby avoid using numerical integration at each step of the iteration.

Here is an outline of that change of variable. You are working in a realm where lambda and c must both be negative with lambda having the greater absolute value. Define A = -c and B = -lambda, so that A and B are positive quantities with A < B. Make the following change of variable:

y = sqrt(B-A)*sin(t)
dy = sqrt(B-A)*cos(t)*dt

(After all the smoke clears away) the first of your integrals becomes the integral with respect to t from 0 to pi/2 of

A/sqrt(A+B-(B-A)*sin(t)^2)

which can be evaluated in terms of a complete elliptic integral of the first kind. It is easy to show that your second integral will reduce to a combination of complete elliptic integrals of the first and second kinds.

Hence in computing your objective function for 'fsolve' you need not call on Matlab's 'integral' function at all but only on the complete elliptic integral functions.