assume can't constrain vpasolve

Hello,

I am attempting to solve a system of nonlinear equations, but unfortunately, vpasolve find answer outside the range I specified and also ignores the constrain condition using "assume".

Here is my code:

% ------------ Parameters --------------------
M = 2;
N = 4;
a = 22.3702/25.4;
beta_10 = 3.41;
k0 = 4.9348;
K2 = 493.2753;
V = [0.7509    1.0000    1.0000    0.7509];
Zb =  [-13.0203-6.9290i  -19.0747-10.8648i  -38.1365-9.9235i  -46.2847-14.8522i];
X =   [0.1672   -0.2434    0.1502   -0.0871 ; 0.0871   -0.1502    0.2434   -0.1672];

% ------------ Required Functions --------------------
syms g(x) h1(y) h2(y) v(x)
g(x) = 1.177.*sin((pi.*x)./a).^2;
h1(y) = 1-275.*(y-1).^2;
h2(y) = -14.*(y-1);
h(y) = h1(y) + 1i*h2(y);
v(x) = 1.517+1.833.*x.^2;

% ------------ Equation Functions --------------------
x = sym('x',[M N]); 
y = sym('y',[M N]);


E1 = []; E2 = []; E3 = []; eq3 = 0; eq4 = 0;
j = 1;
for i = 1:N
    Fn    = ((cos((beta_10/k0)*y(j,i)*v(x(j,i)))-cos(y(j,i)*v(x(j,i))))/sin(y(j,i)*v(x(j,i))))*(sin(pi*x(j,i))/a);
    Ya_G0 = (K2*Fn^2) / (((K2*Fn^2)/(g(x(j,i))*h(y(j,i))))+(Zb(j,i)));

    E1 = [E1; imag((K2*Fn^2)/(g(x(j,i))*h(y(j,i)))) == -1*imag(Zb(j,i))];

    if i==1  && j==1  
         Fn_1   = Fn;
        Ya_G0_1 = Ya_G0;
    else
        
        if mod(i-1,2) == 1
            sign = -1;
        else
            sign = +1;
        end

    E2 = [E2; (Ya_G0/(Fn)) ==  (-1)^(j-1)*sign * abs(V(j,i)/V(1,1)) * (sin(y(j,i)*v(x(j,i)))/(sin(y(1,1)*v(x(1,1))))) * (Ya_G0_1/(Fn_1))   ];
    
    end
    
    eq3 = eq3 + Ya_G0;

end

E3 = [E3; eq3 == 2];


E = [E1.',E2.',E3.'];




assume(x,"real");  assumeAlso(x<0.48*a);
assume(y,"positive"); assumeAlso(y,"real");


X_T = X.';
X_trial = X_T(:).';
sol_initial = [X_trial(1:M*N/2),ones(1,M*N/2)];


% sol_range = [all_x_range(1:M*N/2,:);repmat([0.96 1.04],M*N/2,1)];
% sol_initial = sol_range;

x_T = x.'; x_unknown_vector = x_T(:).';
y_T = y.'; y_unknown_vector = y_T(:).';
unknown = [x_unknown_vector(1:M*N/2),y_unknown_vector(1:M*N/2)];


S = vpasolve(E,unknown,sol_initial);

MyFieldNames = fieldnames(S);


MyValuesX = zeros(1,M*N/2); MyValuesY = zeros(1,M*N/2);
for i=1:M*N/2
MyValuesX(1,i) = getfield(S,MyFieldNames{i});
MyValuesY(1,i) = getfield(S,MyFieldNames{i+N/2});
end

round(MyValuesY,3)

ans = 1×4

0.1500 -0.0870 1.0050 1.0070

I also used the range for "vpasolver" as commented "sol_range" in the above code, but I didn't get the desired result either.

The solution to "y" should be a number around 1.

Thanks for your help in advance.

Answer 1

[Torsten]

vpasolve is a numerical solver - symbolic commands like "assume" are not respected. But you can set ranges for the variables in "vpasolve" by setting "init_param" to a numerical search interval.

[Faraz j]

I defined the range as

all_x_range = repmat([0 a/2;-a/2 0],M*N/2,1);
sol_range = [all_x_range(1:M*N/2,:);repmat([0.96 1.04],M*N/2,1)];

but still give me a solution outside of this range !!!!

[Torsten]

Most probably, no solution for y in the specified range can be found.

[Faraz j]

This is the solution vpasolve should find !

MyValuesX = [0.167 -0.2434 0.1502 -0.0871];
MyValuesY = [1.005  1.007  1.0075  1.0114];

[Torsten]

Ok, vpasolve failed.

Here is a numerical solution with "lsqnonlin":

% ------------ Parameters --------------------
M = 2;
N = 4;
a = 22.3702/25.4;
beta_10 = 3.41;
k0 = 4.9348;
K2 = 493.2753;
V = [0.7509    1.0000    1.0000    0.7509];
Zb =  [-13.0203-6.9290i  -19.0747-10.8648i  -38.1365-9.9235i  -46.2847-14.8522i];
X =   [0.1672   -0.2434    0.1502   -0.0871 ; 0.0871   -0.1502    0.2434   -0.1672];
% ------------ Required Functions --------------------
syms g(x) h1(y) h2(y) v(x)
g(x) = 1.177.*sin((pi.*x)./a).^2;
h1(y) = 1-275.*(y-1).^2;
h2(y) = -14.*(y-1);
h(y) = h1(y) + 1i*h2(y);
v(x) = 1.517+1.833.*x.^2;
% ------------ Equation Functions --------------------
x = sym('x',[M N]); 
y = sym('y',[M N]);
E1 = []; E2 = []; E3 = []; eq3 = 0; eq4 = 0;
j = 1;
for i = 1:N
    Fn    = ((cos((beta_10/k0)*y(j,i)*v(x(j,i)))-cos(y(j,i)*v(x(j,i))))/sin(y(j,i)*v(x(j,i))))*(sin(pi*x(j,i))/a);
    Ya_G0 = (K2*Fn^2) / (((K2*Fn^2)/(g(x(j,i))*h(y(j,i))))+(Zb(j,i)));
    E1 = [E1; imag((K2*Fn^2)/(g(x(j,i))*h(y(j,i)))) - (-1*imag(Zb(j,i)))];
    if i==1  && j==1  
         Fn_1   = Fn;
        Ya_G0_1 = Ya_G0;
    else
        
        if mod(i-1,2) == 1
            sign = -1;
        else
            sign = +1;
        end
        tmp=(Ya_G0/(Fn)) -(  (-1)^(j-1)*sign * abs(V(j,i)/V(1,1)) * (sin(y(j,i)*v(x(j,i)))/(sin(y(1,1)*v(x(1,1))))) * (Ya_G0_1/(Fn_1)) );
    E2 = [E2; tmp  ];
    
    end
    
    eq3 = eq3 + Ya_G0;
end
E3 = [E3; eq3 - 2];
E = [E1.',E2.',E3.'];
Enum = matlabFunction(E);
Enum = @(z)real(Enum(z(1),z(2),z(3),z(4),z(5),z(6),z(7),z(8)));

all_x_range = repmat([0 a/2;-a/2 0],M*N/2,1);
 sol_range = [all_x_range(1:M*N/2,:);repmat([0.96 1.04],M*N/2,1)];
 lb = sol_range(:,1);
 ub = sol_range(:,2);
 format long
 solnum = lsqnonlin(Enum,(lb+ub)/2,lb,ub,optimset('TolX',1e-12,'TolFun',1e-12))

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.

solnum = 8×1

0.167210839303176 -0.243345930008853 0.150185145920223 -0.087101777478330 1.005182169079357 1.007061443401629 1.007465052260141 1.011445713909488

 Enum(solnum)

ans = 1×8

1.0e-12 * -0.071942451995710 0.156319401867222 0.099475983006414 -0.110134124042816 0 -0.000888178419700 -0.001554312234475 0.000444089209850

[Faraz j]

Thanks a lot ! you're life saver

One more issue @Torsten! for different M and N, how can I modify this

Enum = @(z)real(Enum(z(1),z(2),z(3),z(4),z(5),z(6),z(7),z(8)));

[Torsten]

Enum = matlabFunction(E,'Vars',{[x(1,1:N),y(1,1:N)]})
Enum = @(z)real(Enum(z.'));

instead of

Enum = matlabFunction(E);
Enum = @(z)real(Enum(z(1),z(2),z(3),z(4),z(5),z(6),z(7),z(8)));

This is at least correct for the equations you solved above since j was not increased and remained 1 (so I don't know why you created x and y as x(M,N) and y(M,N) and not as x(1,N) and y(1,N))

Answer 2

[Matt J]

There doesn't appear to be any good reason for you to be using sym variables. Your code doesn't seem to make use of any other Symbolic Math Toolbox functions like diff or simplify.... Just re-implement using non-symbolic variables and solve with lsqnonlin which does allow you to impose bounds on the variables. Or, use matlabFunction to convert your symbolic functions into nonsymbolic ones.

[Faraz j]

Thanks for your suggestion !

I wanted to try this approach, but building the function was complicated for me without knowing matlabFunction.

I am not good at this, I think. I tried this and got an error

fun = matlabFunction(E);
lsqnonlin(fun,sol_initial)

what am I doing wrong?

[Matt J]

what am I doing wrong?

Who knows? You haven't shown copy/pastes of the errors.

[Faraz j]

Thanks for your advice, @Torsten solved the issue !