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

조회 수: 2(최근 30일)
Raik Elster 2021년 10월 5일
댓글: Raik Elster 2021년 11월 18일
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
##### 댓글 수: 6표시숨기기 이전 댓글 수: 5
Raik Elster 2021년 10월 6일
편집: Raik Elster 2021년 10월 6일
Now I get it, thank you for your comment. :)
I'll continue the discussion below Björn's answer.

댓글을 달려면 로그인하십시오.

### 답변(2개)

Bjorn Gustavsson 2021년 10월 6일
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
##### 댓글 수: 4표시숨기기 이전 댓글 수: 3
Raik Elster 2021년 10월 7일
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 2021년 10월 7일
편집: Matt J 2021년 10월 7일
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.
##### 댓글 수: 9표시숨기기 이전 댓글 수: 8
Raik Elster 2021년 11월 18일
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!