# Using fsolve to solve equation system for coefficients of a Schwartz-Christoffel map

7 views (last 30 days)
Raik Elster on 5 Oct 2021
Commented: Raik Elster on 18 Nov 2021 at 6:27
Dear MATLAB community,
I'm trying to solve an equation system to get the coefficients of a Schwartz-Christoffel map. In order to understand my problem, I'm giving you a short outline of the underlying mathematics at first.
My goal is to implement a Schwartz-Christoffel map from the z plane to the w plane. The following picture shall give you an impression. The derivation for this map is as follows: The following integrals define then the equation system to be solved for the coefficients A, B, C, D:    Now back to MATLAB. The following code implements the equation system as well as its solution.
clear
% Definition of the geometry
w_slot = 0.120; % Width of the slot distance between rails
h_wg = 0.25; % Height of the electrode rails
w_rail = 0.240; % Width of electrode rails
% Definition of starting point based on coplanar strip waveguide theory
A_0 = abs(w_slot/2 - h_wg/pi); % Using the absolute value here due to the negative sign violating the symmetry
B_0 = w_slot/2;
C_0 = w_slot/2 + w_rail;
D_0 = w_slot/2 + w_rail + h_wg/pi;
x_0 = [A_0;B_0;C_0;D_0];
% Options for the solver
options_fsolve = optimoptions('fsolve','FiniteDifferenceType','central','FunctionTolerance',1e-9,'StepTolerance',1e-12,'OptimalityTolerance',1e-3,'Algorithm','levenberg-marquardt','Display','iter','MaxFunctionEvaluations',1500,'TolPCG',0.001,'MaxPCGIter',6);
% Setup of fsolve
[sol, fval,exflag] = fsolve(@(x) eq_sys(x, w_slot, h_wg, w_rail),x_0,options_fsolve);
% Definition of the derivation of the Schwartz-Christoffel mapping
dz_dw = @(w, x) sqrt((w.^2-x(2).^2).*(w.^2-x(3).^2)./(w.^2-x(1).^2)./(w.^2-x(4).^2));
% Calculation of the integrals in order to verify the found solution
w_slot_test = integral(@(w) dz_dw(w,sol),-sol(1),sol(1));
h_wg_1_test = 2*integral(@(w) dz_dw(w,sol),sol(1),sol(2))/1i;
w_rail_test = integral(@(w) dz_dw(w,sol),sol(2),sol(3));
h_wg_2_test = 2*integral(@(w) dz_dw(w,sol),sol(3),sol(4))/1i;
% Defintion of the equation system to be solved
function S = eq_sys(x, w_slot, h_wg, w_rail)
% Definition of the derivation of the Schwartz-Christoffel mapping
dz_dw = @(w) sqrt((w.^2-x(2).^2).*(w.^2-x(3).^2)./(w.^2-x(1).^2)./(w.^2-x(4).^2));
% Definition of the equation system
S(1) = integral(dz_dw,-x(1),x(1)) - w_slot;
S(2) = integral(dz_dw,x(1),x(2)) - 1i*h_wg/2;
S(3) = integral(dz_dw,x(2),x(3)) - w_rail;
S(4) = integral(dz_dw,x(3),x(4)) + 1i*h_wg/2;
% Definition of the constraints that the solution shall be real due to
% the Schwartz-Christoffel map
S(5) = imag(x(1));
S(6) = imag(x(2));
S(7) = imag(x(3));
S(8) = imag(x(4));
end
The solver stops with exit flag "-2" regardless of the chosen solver (trust-region, trust-region dogleg or levenberg-marquardt). Of course, the integrals to verify the solution don't give the expected values given from the equation system. I also played around with the different tolerance values, which couldn't make the solver work properly. Using lsqnonlin() also didn't work out. Do you have any suggestion about the setup of fsolve or using another solver for this kind of problem? Maybe there could be also a better expression for the equation system. Any kind of help is highly appreciated.
Kind regards,
Raik
Raik Elster on 6 Oct 2021
Now I get it, thank you for your comment. :)
I'll continue the discussion below Björn's answer.

Bjorn Gustavsson on 6 Oct 2021
You can use the optimization-functions to get a solution, or a hopefully good starting-point for fsolve. If you define a sum-of-squares version of your eq_sys like this:
% Defintion of the equation system to be solved
function sumsqrS = sumsqr_eq_sys(x, w_slot, h_wg, w_rail)
% Definition of the derivation of the Schwartz-Christoffel mapping
dz_dw = @(w) sqrt((w.^2-x(2).^2).*(w.^2-x(3).^2)./(w.^2-x(1).^2)./(w.^2-x(4).^2));
% Definition of the equation system
S(1) = integral(dz_dw,-x(1),x(1)) - w_slot;
S(2) = integral(dz_dw,x(1),x(2)) - 1i*h_wg/2;
S(3) = integral(dz_dw,x(2),x(3)) - w_rail;
S(4) = integral(dz_dw,x(3),x(4)) + 1i*h_wg/2;
% Definition of the constraints that the solution shall be real due to
% the Schwartz-Christoffel map
S(5) = imag(x(1));
S(6) = imag(x(2));
S(7) = imag(x(3));
S(8) = imag(x(4));
sumsqrS = sum(real(S(1:4)).^2) + sum(imag(S(1:4)).^2) + sum(S(5:8).^2);
end
Now you can call fminsearch with this function to find an optimum point - and if there are a real x that satisfy all four integrals S(1:4) will all be zero which is a minimum. If no such x is found then you might have gotten a good (better than randomly picked) starting-poin for fsolve to start from. You will definitely get something.
HTH
Raik Elster on 7 Oct 2021
Considering the constraint with fmincon() as well as with minimize() gives a worse solution. Both give the same solution as with fminsearch() when the constraint is omitted. I'll look through the optimization toolbox, if there's another helpful function. I'll also check on Matt's answer below.

Matt J on 7 Oct 2021
Edited: Matt J on 7 Oct 2021
You cannot use complex-valued equations with fsolve, except under specific conditions, so your equations should really be set up this way:
function S = eq_sys(x, w_slot, h_wg, w_rail)
% Definition of the derivation of the Schwartz-Christoffel mapping
dz_dw = @(w) sqrt((w.^2-x(2).^2).*(w.^2-x(3).^2)./(w.^2-x(1).^2)./(w.^2-x(4).^2));
% Definition of the equation system
T(1) = integral(dz_dw,-x(1),x(1)) - w_slot;
T(2) = integral(dz_dw,x(1),x(2)) - 1i*h_wg/2;
T(3) = integral(dz_dw,x(2),x(3)) - w_rail;
T(4) = integral(dz_dw,x(3),x(4)) + 1i*h_wg/2;
S(1:4)=real(T);
S(5:8)=imag(T);
end
Note however, that your integrals are improper because some of the limits of integration lie at the poles of dz/dw, so numerical integration via integral() is going to be unrelieable.
Raik Elster on 18 Nov 2021 at 6:27
Finally, after some hours of thinking I got a working solution. The main problem was the definition of dz_dw. At some point using the sqrt() function affects the calculation of the integral at solving the equation system as well as at verifying the solution. By replacing sqrt() with the exponent 0.5 or -0.5 respectively the solution becomes consistent and is verified by the integrals using the found solution. In the end, the following expression has to be used:
dz_dw = @(w) (w.^2-x(2).^2).^0.5.*(w.^2-x(3).^2).^0.5.*(w.^2-x(1).^2).^(-0.5).*(w.^2-x(4).^2).^(-0.5);

R2020a

### Community Treasure Hunt

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

Start Hunting!