2 views (last 30 days)

Show older comments

Faraz j on 4 Feb 2024

Edited: Torsten on 6 Feb 2024

Accepted Answer: Torsten

Open in MATLAB Online

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.

##### 0 Comments Show -2 older commentsHide -2 older comments

Show -2 older commentsHide -2 older comments

Sign in to comment.

Sign in to answer this question.

### Accepted Answer

Moved: Torsten on 4 Feb 2024

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.

##### 6 Comments Show 4 older commentsHide 4 older comments

Show 4 older commentsHide 4 older comments

Faraz j on 4 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3053106

Open in MATLAB Online

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 on 4 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3053306

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

Faraz j on 5 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3053696

Edited: Faraz j on 5 Feb 2024

Open in MATLAB Online

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 on 5 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3053791

Edited: Torsten on 5 Feb 2024

Open in MATLAB Online

Ok, vpasolve failed.

Here is a numerical solution with "lsqnonlin":

% ------------ Parameters --------------------

M = 2;

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 on 6 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3054861

Edited: Faraz j on 6 Feb 2024

Open in MATLAB Online

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 on 6 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3055721

Edited: Torsten on 6 Feb 2024

Open in MATLAB Online

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

Sign in to comment.

### More Answers (1)

Matt J on 4 Feb 2024

Edited: Matt J on 4 Feb 2024

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.

##### 3 Comments Show 1 older commentHide 1 older comment

Show 1 older commentHide 1 older comment

Faraz j on 5 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3053721

Edited: Faraz j on 5 Feb 2024

Open in MATLAB Online

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 on 5 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3053941

what am I doing wrong?

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

Faraz j on 6 Feb 2024

#### Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/2078046-assume-can-t-constrain-vpasolve#comment_3054866

Thanks for your advice, @Torsten solved the issue !

Sign in to comment.

Sign in to answer this question.

### See Also

### Categories

Mathematics and OptimizationSymbolic Math ToolboxSymbolic Computations in MATLABConversion Between Symbolic and Numeric

Find more on **Conversion Between Symbolic and Numeric** in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français

- United Kingdom(English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)

Contact your local office