MATLAB Answers

## Solving a system of integral equations numerically

Asked by Lewis Hancox

### Lewis Hancox (view profile)

on 12 Mar 2018
Latest activity Commented on by Star Strider

### Star Strider (view profile)

on 12 Mar 2018
Accepted Answer by Star Strider

### Star Strider (view profile)

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

#### 0 Comments

Sign in to comment.

## 3 Answers ### Star Strider (view profile)

Answer by Star Strider

### Star Strider (view profile)

on 12 Mar 2018
Accepted Answer

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.

Star Strider

### Star Strider (view profile)

on 12 Mar 2018
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

### Lewis Hancox (view profile)

on 12 Mar 2018
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

### Star Strider (view profile)

on 12 Mar 2018
As always, my pleasure!

Sign in to comment. ### Roger Stafford (view profile)

Answer by Roger Stafford

### Roger Stafford (view profile)

on 12 Mar 2018

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.

#### 0 Comments

Sign in to comment. ### Roger Stafford (view profile)

Answer by Roger Stafford

### Roger Stafford (view profile)

on 12 Mar 2018
Edited by Roger Stafford

### Roger Stafford (view profile)

on 12 Mar 2018

@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.

#### 0 Comments

Sign in to comment.