How to solve a system of equations whose unknowns values are exponents with cole-cole methods

I want to solve this equation,

I want to calculate the value of B1,B2

I use this equation

eh = 45.;
kl = 4.15;
de1 = 6200;
fc1 = 500*10^3;
de2 = 990;
fc2 = 2.70*10^6;
b1 = 0.86;
b2=0.98;
f=10^4;
% ev=8.854187817 * 10^(-12);
ev=8.854 * 10^(-12); %pF/m

e0 = eh+de1/(1+j*f/fc1)^b1+de2/(1+j*f/fc2)^b2+kl/(j*2*pi*f*ev)
e1 = real(e0);
e2 = -imag(e0);

%%
x0 = [b1,b2];
% options=optimset('MaxFunEvals',1e4,'MaxIter',1e4);
options = optimoptions('fsolve','Display','iter','TolFun',1e-30,'TolX',1e-30);

[x,FVAL,EXITFLAG,OUTPUT,JACOB] = fsolve(@(x)fun_cole_cole(x,fc1,fc2,f,de1,de2,eh,kl,ev,e1,e2),x0,options);

The fun_cole_cole function is shown below

function f = fun_cole_cole(x,fc1,fc2,f,de1,de2,eh,kl,ev,e1,e2)  % b c可以是随意的参数



A1 = fc1^x(1)+f^x(1)*cos(pi*x(1)/2);
A2 = fc2^x(2)+f^x(2)*cos(pi*x(2)/2);
B1 = f^x(1)+sin(pi*x(1)/2);
B2 = f^x(2)+sin(pi*x(2)/2);
D1 = de1*fc1^x(1);
D2 = de2*fc2^x(2);

w = 2*pi*f;

f1 = A1*D1/(A1^2+B1^2)+A2*D2/(A2^2+B2^2)+eh-e1;
f2 = B1*D1/(A1^2+B1^2)+B2*D2/(A2^2+B2^2)+kl/(w*ev)-e2;

f = [f1;f2];

But I got this result

No solution found.

fsolve stopped because the problem appears to be locally singular.

When I set options as

options=optimset('MaxFunEvals',1e4,'MaxIter',1e4);

I got this

No solution found.

fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the default value of the function tolerance.

Answer 1

[Walter Roberson]

This reacts well to a standard trick of converting a system of equations into a least-squared minimization.

The fsolve() is not needed; I kept it in to show you that the result of the fmincon() checks out with fmincon()

format long g
eh = 45.;
kl = 4.15;
de1 = 6200;
fc1 = 500*10^3;
de2 = 990;
fc2 = 2.70*10^6;
b1 = 0.86;
b2=0.98;
f=10^4;
% ev=8.854187817 * 10^(-12);
ev=8.854 * 10^(-12); %pF/m

e0 = eh+de1/(1+j*f/fc1)^b1+de2/(1+j*f/fc2)^b2+kl/(j*2*pi*f*ev)

e0 = 
           7233.00405024904 -      7459936.41901046i

e1 = real(e0);
e2 = -imag(e0);

%%
x0 = [b1,b2];
% options=optimset('MaxFunEvals',1e4,'MaxIter',1e4);
options = optimoptions('fsolve','Display','iter','TolFun',1e-3,'TolX',1e-3);
f = @(x)fun_cole_cole(x,fc1,fc2,f,de1,de2,eh,kl,ev,e1,e2);

syms x [1 2]

residue = sum(f(x).^2);
R  = matlabFunction(residue,'vars',{x});

[best, fval] = fmincon(R, x0)

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.

best = 1×2

1.09674543072127 0.653154719676118

fval = 
      1.68411061989985e-11

[x,FVAL,EXITFLAG,OUTPUT,JACOB] = fsolve(f,best,options)

                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3     1.68383e-11                       0.00096               1

Equation solved at initial point.

fsolve completed because the vector of function values at the initial point
is near zero as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

x = 1×2

1.09674543072127 0.653154719676118

FVAL = 2×1

2.15581258089514e-06 3.49152833223343e-06

EXITFLAG = 
     1

OUTPUT = struct with fields:
       iterations: 0
        funcCount: 3
        algorithm: 'trust-region-dogleg'
    firstorderopt: 0.000960191844270147
          message: '↵Equation solved at initial point.↵↵fsolve completed because the vector of function values at the initial point↵is near zero as measured by the value of the function tolerance, and↵the problem appears regular as measured by the gradient.↵↵<stopping criteria details>↵↵Equation solved. The final point is the initial point.↵The sum of squared function values, r = 1.683830e-11, is less than sqrt(options.FunctionTolerance) = 3.162278e-02.↵The relative norm of the gradient of r, 9.601918e-04, is less than options.OptimalityTolerance = 1.000000e-03.↵↵'

JACOB = 2×2

91.0225376370537 112.5927734375 -331.207215298002 -134.375

function f = fun_cole_cole(x,fc1,fc2,f,de1,de2,eh,kl,ev,e1,e2)  % b c可以是随意的参数

A1 = fc1^x(1)+f^x(1)*cos(pi*x(1)/2);
A2 = fc2^x(2)+f^x(2)*cos(pi*x(2)/2);
B1 = f^x(1)+sin(pi*x(1)/2);
B2 = f^x(2)+sin(pi*x(2)/2);
D1 = de1*fc1^x(1);
D2 = de2*fc2^x(2);

w = 2*pi*f;

f1 = A1*D1/(A1^2+B1^2)+A2*D2/(A2^2+B2^2)+eh-e1;
f2 = B1*D1/(A1^2+B1^2)+B2*D2/(A2^2+B2^2)+kl/(w*ev)-e2;

f = [f1;f2];
end

[HONG CHENG]

Dear sir, thank you for your answer.

My matlab is R2018a, when I run your codes, I have errors,

Error using syms (line 222)
Invalid variable name.

Error in cole_cole_Test2 (line 31)
syms x [1 2]

Here are all the codes

format long g

eh = 45.;
kl = 4.15;
de1 = 6200;
fc1 = 500*10^3;
de2 = 990;
fc2 = 2.70*10^6;
b1 = 0.86;
b2=0.98;
f=10^4;
% ev=8.854187817 * 10^(-12);
ev=8.854 * 10^(-12); %pF/m

e0 = eh+de1/(1+j*f/fc1)^b1+de2/(1+j*f/fc2)^b2+kl/(j*2*pi*f*ev)


%%

e1 = real(e0);
e2 = -imag(e0);

x0 = [b1,b2];
options = optimoptions('fsolve','Display','iter','TolFun',1e-3,'TolX',1e-3);

f = @(x)fun_cole_cole(x,fc1,fc2,f,de1,de2,eh,kl,ev,e1,e2);

syms x [1 2]
residue = sum(f(x).^2);
R  = matlabFunction(residue,'vars',{x});
[best, fval] = fmincon(R, x0)

This line, I am not familar with this mode

syms x [1 2]

Is this a new feture of new versions of Matlab?

[Walter Roberson]

Replace

syms x [1 2]

with

x = sym('x', [1 2])

[HONG CHENG]

Thank you for your helps

[HONG CHENG]

I think I have another question. I am not sure if my understanding is right or not.

we can use

[best, fval] = fmincon(R, x0)

to finds a constrained minimum of a function of several variables.

after we get the values of best, we use them as the intial values for fsolve

options = optimoptions('fsolve','Display','iter','TolFun',1e-3,'TolX',1e-3);
[x,FVAL,EXITFLAG,OUTPUT,JACOB] = fsolve(f,best,options)

But finally, we get the same results and also get the values of EXITFLAG,OUTPUT,JACOB

[Walter Roberson]

As I wrote, "The fsolve() is not needed; I kept it in to show you that the result of the fmincon() checks out with fmincon()"

fmincon() is a minimizer; it could potentially have returned a value that was not a root, or it might have gotten stuck just a little away from a root. Any time you switch between techniques (I recoded fsolve as fmincon) then at the beginning you want to reassure yourself that the new technique does as well as the old one would have done. Once you have tried a few cases and seen that it does work, you can get rid of the fsolve call.