Solving non-linear integral system with complex numbers

6 views (last 30 days)
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!
  2 Comments
Eggstrema
Eggstrema on 11 Oct 2021
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.

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 10 Oct 2021
Edited: Walter Roberson on 10 Oct 2021
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"
  4 Comments
Walter Roberson
Walter Roberson on 11 Oct 2021
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"

Sign in to comment.

More Answers (0)

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!