Using fsolve to solve equation system for coefficients of a Schwartz-Christoffel map
1 view (last 30 days)
Show older comments
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 Comments
nathan blanc
on 6 Oct 2021
Raik, Like bjorn said: if you have four functions (A,B,C,D) that you want to minimize define your cost function as A^2+B^2+C^2+D^2, and you can minimize all of them together.
Answers (2)
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
4 Comments
Bjorn Gustavsson
on 7 Oct 2021
When I need to solve optimization-problems with constraints on the parameters I use fminsearchbnd from the file exchange (there is also minimize which might be a more powerful version), if you have access to fmincon then that is the built-in matlab-function for constrained optimization. Then you can force A B C and D to be in accending order.
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.
9 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!