Roots of a two variables equation involving a numerical integration

I'm looking for a solution of a of two variables equation problem,

The issue is that '' and '' are integrals of a certain kind of functions 'f' and 'g' involving the variables x, y and p, in a very complicate way (since the integration is only on p, it eventually disappears), that include logarithms of the three variables, like the example function below

I have tried to use a double for loop for the variables x and y, so that I can integrate numerically over p only. However, it has not been efficient and one of the problems I have is that for certain configurations of values of x,y and p, negative numbers are obtained that generate an indeterminacy in the logarithm.

In summary, I must find the pairs of numbers (x, y) that satisfy the following type of equation, and plot them

[Ameer Hamza]

Can you show us the functions f and g? You can try fsolve(), or even fmincon() if there are some bounds on the values of x and y.

[Daniel Tamayo]

Thanks for your answer! sorry about the length...

Where,

and

and

[Ameer Hamza]

Yes, these are quite long. If you can write them in MATLAB, maybe we can suggest a solution with fsolve, fmincon or some other optimization algorithm.

Answer 1

[J. Alex Lee]

Assuming x and y are scalars and there is a solution, use fsolve, not an optimization algorithm. You'll need initial guesses. The beauty of numerical solution is that you can literally just define the functions you have written out, and call them. The numerical integration can be done with matlab's in-built integral() function, which according to docs can even take Inf as an integration limit (sounds cool, i have never tried this). The structure would look something like

X0 = [x0,y0]

X = fsolve(@(X)myresidual(X),X0)

function res = myresidual(X)
x = X(1);
y = X(2);
res = (y-x) .* (1+integral(@(p)f(x,y,p),0,Inf)) + integral(@(p)g(x,y,p),0,Inf);
end

function out = f(x,y,p)
% and you call other functions to do the other complicated sub-parts like eps1, L1p, L1m, etc.
end
function out = g(x,y,p)

end

function out = L1(x,y,M,e1)
end

% etc

[Daniel Tamayo]

Thanks for your help J. Alex Lee!!

But, the problem I see is that ‘fsolve’ solves a nxn system, so ‘fsolve’ only gives one root.

However, I need all the points (x_i,y_i) that are solutions of the equation (in a finite range of x and y, because there are infinite solutions), and plot them.

Do you know if there is any way to find all the roots of a system in which there is an equation with two unknowns?

Thanks again!

[J. Alex Lee]

The fact that fsolve solves a system of nxn equations has nothing to do with the fact that it finds only 1 root, unless I'm misunderstanding your statement.

In any case, for nonlinear systems, there is no way (in general) to even know how many roots there are (or if there are any). As far as I know there are no general numerical methods that attempt to "find all" roots.

Since you seem to know there are infinitely many solutions, do you roughly where/how they are spaced? If you know roughly where they are, you can provide them as initial guesses to fsolve so that it will converge to that root.If you know spacing, you can find 1 root, then adjust it by the known spacing and re-solve, and so on.

[John D'Errico]

I think the problem you are tripping over is that fsolve returns only one solution. But here, you have apparently two variables and one equation. (Your commented response to Ameer is a bit confusing though.)

That means there will generally be infinitely many solutions, think of it as a curve in the (x,y) plane. This is often called an implicit function, much like a circle in the plane.

But finding all solutions if an analytical solution fails to exist is impossible, since there are infinitely many of them. fsolve can never do that, since it can solve for only one solution at a time, based on the starting values to the problem.

If you just want a plot, then fimplicit is a great way to solve these problems.

IF you just want a list of a FINITE set of solutions for whatever reason, then one thing is to just write a loop, fixing x at a sequence of values. Then solve for y, as essentially a function of x. Save all solutions found for y as a function of x.

A problem may arise with this approach, even for a simple case - the circle. Here you would find for any given x, there are 0, 1 or 2 solutions for y. The real world sucks, and sometimes there is no simple answer to a problem.

Another solution might be to use a contour plot, here the proper tool will be contourc. The idea of contourc is it will have you evaluate the overall function for a large grid of points. Then you will tell it to search for the piecewise linear curve through that grid where the function is exactly zero (determined by a local linear interpolation.) The result will be a segmented linear approximation to the solutions of your problem, in the (x,y) plane. The quality of this approimation will be a function of how fine was your grid as provided to contourc.

[J. Alex Lee]

Yes, my oversight, thanks John D'Errico, I agree with solutions proposed, contour being an interesting solution.

To the solution of fixing x and finding y, you can use path following strategies that are immune to turning points, such as [pseudo] arc-length continuation. There are more advanced methods for dealing with other types of singularties like bifurcations, but i'm not as familiar.

[John D'Errico]

One thing that I love is the ability of contour and contourc to solve problems of this sort, as well as return the solution locus as an approximation.

Path following is another nice idea, one that is also often overlooked. I even considered writing a toolbox for the purpose once.

[J. Alex Lee]

Indeed, I am using contourX more and more in lieu of actually solving for conditions when the problem is challenging. ironically i am hell bent on analytical solutions when my problem reduces to a cubic polynomial.

I never bothered to google it until now, but there's a nice wiki page with links to toolboxes in different languages

https://en.wikipedia.org/wiki/Numerical_continuation

including 1 for matlab, though i haven't at all looked at it:

https://sourceforge.net/projects/matcont/