Solving non-linear integral system with complex numbers

I am trying to solve this nonlinear system of equations numerically:

with the following code:

% param
W = 1e-3;
S = 0.5e-3;
t = 0.1e-3;

% initial guess
a0 = W/2;
b0 = W/2+t/pi;
c0 = W/2+S-t/pi;
d0 = W/2+S;
s0 = [a0,b0,c0,d0];

% eqs and solve
fun1 = @(s) integral(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),-s(1),s(1))-W;
fun2 = @(s) integral(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(1),s(2))+1i*t/2;
fun3 = @(s) integral(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(2),s(3))-S;
fun4 = @(s) integral(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(3),s(4))-1i*t/2;
s = fsolve(@(s) [fun1(s),fun2(s),fun3(s),fun4(s)],s0)

However, I am getting unexpected results because the relationship of the solutions should be and that they are all real-valued. The command line also tells me that that "equation solved, solver stalled" or "equation solved, possible inaccuracies" depending on my input parameters W, S, and t. Any help would be appreciated, thanks!

[Walter Roberson]

I just realized that you do not appear to have stated that the s values are strictly positive. You said they are real-valued and s1 < s2 < s3 < s4, but is it possible for them to be negative?

[Eggstrema]

This is based on a Schwarz–Christoffel conformal mapping problem, where I am mapping points mostly to the first quadrant so I am expecting positive values. It is possible for the solutions to be negative (second quadrant), but the constraints would also change to s4 < s3 < s2 < s1 since it will be reflected on the imaginary axis.

Answer 1

[Walter Roberson]

fsolve() is completely unable to handle constraints.

format long g
% param
W = 1e-3;
S = 0.5e-3;
t = 0.1e-3;

% initial guess
a0 = W/2;
b0 = W/2+t/pi;
c0 = W/2+S-t/pi;
d0 = W/2+S;
s0 = [a0,b0,c0,d0];

% eqs and solve
fun1 = @(s) integral(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),-s(1),s(1))-W;
fun2 = @(s) integral(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(1),s(2))+1i*t/2;
fun3 = @(s) integral(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(2),s(3))-S;
fun4 = @(s) integral(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(3),s(4))-1i*t/2;
residue = @(s) norm([fun1(s), fun2(s), fun3(s), fun4(s)]);
A = [1 -1 0 0; 0 1 -1 0; 0 0 1 -1]; b = eps*ones(3,1);
Aeq = []; beq = [];
lb = zeros(1,4); ub = [];
[bests, fval] = fmincon(residue, s0, A, b, Aeq, beq, lb, ub);

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

string(bests)

ans = 1×4 string array
    "0.00052471"    "0.00054604"    "0.00090234"    "0.00094363"

string(fval)

ans = "8.8712e-05"

[Eggstrema]

Thanks for your answer, is there a reason why the b vector is eps instead of zero?

[Walter Roberson]

When you use an A b pair, the logic is A*x <= b must be satisfied.

You want s1 < s2 with a strict inequality. It needs to be a strict inequality because of the divisions you do: an equality could result in a division by 0.

s1 < s2 implies s1 - s2 < 0 implies [1 -1 0 0] * [s1;s2;s3;s4] < 0 so that would be an A row of [1 -1 0 0] and a b of 0 ... except for the fact that the logic is <= b but you need < b instead. So we nudge slightly using a small number instead of exact 0.

However, it is worth thinking more carefully about whether it should be +eps or -eps

s1 < s2 implies s1 + delta = s2 for some strictly positive delta, implies s1 - s2 = -delta, again strictly positive delta

So what I should have coded is

b = -eps*ones(3,1);

You could use a different number with absolute value less than eps if you wanted, such as -realmin*ones(3,1)

[Eggstrema]

Thanks for the explanation. I plugged back in the solution into the original four equations to see whether I get back W, S, and +/-t/2.

% verify

I1 = fun1(s);

I2 = fun2(s);

I3 = fun3(s);

I4 = fun4(s);

disp([I1+W I2-j*t/2 I3+S I4+j*t/2]);

Columns 1 through 2

    1.00029318882667e-003 + 0.00000000000000e+000i    0.00000000000000e+000 + 6.46803411013140e-006i

  Columns 3 through 4

    500.940271615251e-006 + 0.00000000000000e+000i    0.00000000000000e+000 + 30.2429802599424e-006i

Are there ways to improve the solution? I tried chaging the algorithm and tolerances but to no avail.

[Walter Roberson]

I realized that the integration parts of fun2 and fun4 produced pure real numbers times sqrt(-1), and then you have a pure imaginary part being added to that. That is equivalent to sqrt(-1) times (sum of two real). When you norm() that, then the sqrt(-1) disappears . instead of trying to integrate with complex parts, you can integrate the negative of the parts (inside the sqrt) and remove the 1i* from the constant, transforming the two parts into pure real, which is easier to deal with.

In the below, I am not certain why I am getting imaginary parts in the back-substitutions, especially since fun3 should be sqrt(positive*negative / (negative * positive)) which should simplify to sqrt(positive) with no possibility for complex parts.

If you set digits to 18 or more, then the vpa() in the back substitution fails, implying that it could not meet precision goals. I do not know at the moment why that would be happening.

When you compare fval here the previous, remember that the previous is sqrt() relative to this one -- the one below is equivalent to norm^2

W = 1e-3;
S = 0.5e-3;
t = 0.1e-3;
% initial guess
a0 = W/2;
b0 = W/2+t/pi;
c0 = W/2+S-t/pi;
d0 = W/2+S;
s0 = [a0,b0,c0,d0];
% eqs and solve
fun1 = @(s) int(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),-s(1),s(1))-W;
fun2 = @(s) int(@(z) sqrt(-(z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(1),s(2))+t/2;
fun3 = @(s) int(@(z) sqrt((z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(2),s(3))-S;
fun4 = @(s) int(@(z) sqrt(-(z.^2-s(1).^2).*(z.^2-s(4).^2)./((z.^2-s(2).^2).*(z.^2-s(3).^2))),s(3),s(4))-t/2;

residue = @(s) fun1(s).^2 + fun2(s).^2 + fun3(s).^2 + fun4(s).^2;
syms S [1 4] real
R = matlabFunction(residue(S), 'vars', {S})

R = function_handle with value:
    @(in1)(integral(@(z)sqrt(((in1(:,1).^2-z.^2).*(in1(:,4).^2-z.^2))./((in1(:,2).^2-z.^2).*(in1(:,3).^2-z.^2))),-in1(:,1),in1(:,1))-1.0./1.0e+3).^2+(integral(@(z)sqrt(((in1(:,1).^2-z.^2).*(in1(:,4).^2-z.^2))./((in1(:,2).^2-z.^2).*(in1(:,3).^2-z.^2))),in1(:,2),in1(:,3))-5.0e-4).^2+(integral(@(z)sqrt(-((in1(:,1).^2-z.^2).*(in1(:,4).^2-z.^2))./((in1(:,2).^2-z.^2).*(in1(:,3).^2-z.^2))),in1(:,1),in1(:,2))+5.0e-5).^2+(integral(@(z)sqrt(-((in1(:,1).^2-z.^2).*(in1(:,4).^2-z.^2))./((in1(:,2).^2-z.^2).*(in1(:,3).^2-z.^2))),in1(:,3),in1(:,4))-5.0e-5).^2

A = [1 -1 0 0; 0 1 -1 0; 0 0 1 -1]; b = -realmin*ones(3,1);
Aeq = []; beq = [];
lb = zeros(1,4); ub = [];
[bests, fval] = fmincon(R, s0, A, b, Aeq, beq, lb, ub)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

bests = 1×4

0.0009 0.0011 0.0015 0.0016

fval = 4.5751e-07

string(bests)

ans = 1×4 string array
    "0.00085322"    "0.0010892"    "0.0014844"    "0.0015651"

digits(17)

vpa(fun1(sym(bests)))

ans = 

0.00022868074369929651

vpa(fun2(sym(bests)))

ans = 

0.00043562490278369904

vpa(fun3(sym(bests)))

ans = 

vpa(fun4(sym(bests)))

ans = 

string(ans)

ans = "0.00010347939488391012 + 0.000000000000000068181053664877063i"