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

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

[nathan blanc]

you could try using the optimization toolbox or global optimization toolbox, where your cost function is the differrence between the left and right sides of your equations.

[Raik Elster]

Thanks for your answer. I had a small glimpse into optimization and my impression was, that it only works for scalar optimization functions and not for multidimensional ones as in this case. Please correct me, if I'm wrong.

[Bjorn Gustavsson]

If you look at the sum of the squares (of the real and imaginary parts) of your S(1), S(2), S(3) and S(4) you will definitely have a minimum for an exact solution - if one exists and the optimization-function manages to find it...

[Raik Elster]

I don't doubt that the optimization function does find a solution for one of the equations S(i). This is how I understand the optimization toolbox; that it works for only one function at a time. But in this equation system all the equations S(i) have to be solved at the same time, for which I couldn't find a suitable approach when using optimization. As I mentioned before, please correct me, if I'm wrong with my impression.

[nathan blanc]

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.

[Raik Elster]

Now I get it, thank you for your comment. :)

I'll continue the discussion below Björn's answer.

Answer 1

[Bjorn Gustavsson]

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

[Bjorn Gustavsson]

...and as always: Have a look on the file exchange to see if any of the contributions for this type of task solves your problem: schwarzchristoffel

[Raik Elster]

Your suggestion improved the solution a lot, thank you very much. However, it's still not completely consistent, i.e. the solution for x(3) is larger than x(4), which presumably leads to the abs(h_wg_2_test) being only the half of h_wg_1_test. Both are supposed to be identical.

I implemented fminsearch with your proposed function as follows including an updated fsolve command:

% Options for the solver

options_fminsearch = optimset('Display','iter','TolFun',1e-6,'TolX',1e-6);

% Setup of fminsearch

[sol_fmin, fval_fmin, exflag_fmin] = fminsearch(@(x) sumsqr_eq_sys(x, w_slot, h_wg, w_rail),x_0,options_fminsearch);

% Setup of fsolve

[sol, fval,exflag] = fsolve(@(x) eq_sys(x, w_slot, h_wg, w_rail),sol_fmin ,options_fsolve);

In the end, the solution remains the same after using fsolve. Do you have an idea, how to improve the way of finding a better solution or at least how to check, if it's the best one? I used also rand(1,4) instead of x_0 in fminsearch and fval_fmin won't go below 0.0078, which belongs to the solution mentioned above.

[Bjorn Gustavsson]

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.

[Raik Elster]

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.

Answer 2

[Matt J]

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]

Thank you for your answer. With changing the function like you suggested fsolve() gives a similar solution like fminsearch(), as Björn proposed above. It seems likely that it could be an issue of the integration at the poles. Do you have an idea how to make the integrals more reliable?

[Matt J]

It also appears that you have a constraint 0<=A<B<C<D that needs to be enforced. You can do that with fmincon or, as below, with lsqnonlin. The lsqnonlin solution requires reparametrizing the problem, such that your unknowns are now dx=[A, B-A, C-B, D-C]>0 and x=cumsum(dx).

This gives a solution that satisfies well all but the last equation

>> residual =

Columns 1 through 7

6.6886e-09 0 -1.1171e-10 0 0 4.4024e-09 0

Column 8

0.125

so you probably need a better initial guess for dx_0(4)

% 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_lsqnonlin =  optimoptions('lsqnonlin','FiniteDifferenceType','central',...
                  'FunctionTolerance',1e-12,'StepTolerance',1e-12,'OptimalityTolerance',...
                  1e-12,'Display','iter','MaxFunctionEvaluations',1e6,'MaxIter',1e3 );
% Setup of lsqnonlin
dx_0=[A_0,B_0, w_rail, h_wg/pi];
fun=@(dx) eq_sys(dx, w_slot, h_wg, w_rail);
[dsol, fval,residual, exflag] = lsqnonlin(fun,dx_0,[0,0,0,0]+1e-8,[],options_lsqnonlin);

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          9       0.0753721                        0.0676
     1         18       0.0294003       0.227818        0.00645      
     2         27        0.017674       0.166091        0.00194      
     3         36       0.0156505      0.0783198       0.000689      
     4         45       0.0156253      0.0096488       6.24e-05      
     5         54        0.015625     0.00110729       2.35e-05      
     6         63        0.015625    1.04912e-05       4.32e-09      
     7         72        0.015625     1.4905e-07       4.32e-09      
     8         81        0.015625    3.72626e-08       4.32e-09      
     9         90        0.015625    9.31564e-09       4.32e-09      
    10         99        0.015625    2.32891e-09       4.32e-09      
    11        108        0.015625    5.82228e-10       4.32e-09      
    12        117        0.015625    1.45557e-10       4.32e-09      
    13        126        0.015625    3.63892e-11       4.32e-09      
    14        135        0.015625    9.09731e-12       3.18e-09      

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

sol=cumsum(dsol);


function S = eq_sys(dx, w_slot, h_wg, w_rail)
    
    x=cumsum(dx);

    % 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

[Raik Elster]

Thank you very much for your suggestion, this looks promising. I'm checking on the initial guess and see, if there's a better one.

At least I found a solution to implement the constraint 0<A<B<C<D with fsolve() by adding the following equations to S:

S(9) = sign(x(4)-x(3)) - 1;

S(10) = sign(x(3)-x(2)) - 1;

S(11) = sign(x(2)-x(1)) - 1;

S(12) = sign(x(1)) - 1;

[Matt J]

At least I found a solution to implement the constraint 0<A<B<C<D with fsolve() by adding the following equations to S:

I really wouldn't recommend that. fsolve and lsqnonlin both assume your equation function is continuously differentiable. Since sign() has a discontinuity in it, this will not be the case in your formulation.

[Bjorn Gustavsson]

@Raik Elster: It is way prefereable to define your parameters as Matt J hinted at above:

pars = [xA,dxAtoB,dxBtoC,dxCtoD]

Then in your function fun you calculate the break-points A, B, C and D:

ABCD = cumsum(pars); % dx in Matt J's example above

Now you can use the constrained optimization-functions with the constraint that all four elements in pars should be positive.

One further point to consider is to perhaps make the integration-limits to not include the points where the denominators equals zero - stop some small eps before that. Some care might have to be taken when choosing that...

[Raik Elster]

@Matt J: Thanks for the remark regarding the sign() function. That sounds reasonable to me.

@Bjorn Gustavsson: I played around a litte bit with the limits of the integrals, which didn't make any improvements.

Actually, I tried to have a deeper look into the equation system and it turned out that the fourth equation has to be changed regarding the sign of the result of the integral in order to give a consistent solution:

That's the log of lsqnonlin():

Norm of First-order

Iteration Func-count f(x) step optimality

0 9 0.0208308 0.0939

1 18 0.00215207 0.155293 0.00387

2 27 0.000353583 0.0620739 0.000775

3 36 2.8857e-05 0.0314374 0.000648

4 45 2.66821e-07 0.0115146 0.000106

5 54 9.22245e-12 0.000738071 2.31e-05

6 63 4.15361e-21 3.32849e-06 5.1e-11

7 72 1.18083e-27 9.3658e-11 4.97e-14

Local minimum found.

Optimization completed because the size of the gradient is less than

the value of the optimality tolerance.

This leads to the following solution, which looks very good:

fval =

1.1808e-27

sol =

0.0026 0.0461 0.4066 0.4903

exflag =

1

Apparently, this is contrary to my understanding of the Schwarz-Christoffel map, since the integral goes along the negative direction of the axis of the imaginary part in the z plane, hence I expect a negative sign in the result. On the other hand, the actual mapping is defined by the derivative, where the angles and directions of the map are already considered. Do you have an explanation for that?

[Matt J]

This leads to the following solution, which looks very good:

Glad it's working, but please Accept-click the answer, since it appears to have addressed your problem.

On the other hand, the actual mapping is defined by the derivative, where the angles and directions of the map are already considered.

You're using hte same integrand in all 4 cases. I don't see how an integrand could ever be aware of what direction it's being integrated along and flip its own sign accordingly...

[Raik Elster]

You're using hte same integrand in all 4 cases. I don't see how an integrand could ever be aware of what direction it's being integrated along and flip its own sign accordingly...

You're right, that's the reason why I don't get that a solution is found by simply changing the sign. When talking about the derivative, I was imagining how the integration curve would look like in the z plane. Based on the new equation, the integration goes upwards to i*h_wg instead of going downwards to 0. The derivative, however, is defined by the angels at the corresponding points according to the definition of the Schwarz-Christoffel map. Due to the discrepancy between the integral curve and the mapping derivative based on the actual curve I suppose that there might be a certain mathematical property I don't know of yet.

[Raik Elster]

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);