Solving 4 nonlinear equation for 4 unknowns numerical - doesn't find an answer [0×1 sym]

Hello,

I am attempting to solve a system of nonlinear equations, but unfortunately, vpasolve couldn't find any solutions.

Here is my code:

% ------------ Parameters --------------------
a = 0.9;
beta_10 = 3.5668;
k0 = 4.9907;
K2 = 480;

% ------------ 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 --------------------
syms x1 x2 y1 y2

Fn1 = ((cos((beta_10/k0)*y1*v(x1))-cos(y1*v(x1)))/sin(y1*v(x1)))*(sin(pi*x1)/a);
Fn2 = ((cos((beta_10/k0)*y2*v(x2))-cos(y2*v(x2)))/sin(y2*v(x2)))*(sin(pi*x2)/a);

Ya1_G0 = (K2*Fn1^2) / (((K2*Fn1^2)/(g(x1)*h(y1)))+(3.0691-13.7872i));
Ya2_G0 = (K2*Fn2^2) / (((K2*Fn2^2)/(g(x2)*h(y2)))+(-1.9145-11.3728i));

% ------------ Equations --------------------
E1 = imag((K2*Fn1^2)/(g(x1)*h(y1))) == 13.7872;
E2 = imag((K2*Fn2^2)/(g(x2)*h(y2))) == 11.3728;
E3 = Ya2_G0/Fn2 == -2 * (Ya1_G0/Fn1);
E4 = Ya1_G0 + Ya2_G0 == 0.5;

E = [E1,E2,E3,E4];
S = vpasolve(E,[x1,x2,y1,y2],[-a/4 a/4;-a/4 a/4;0.98 1.02;0.98 1.02])

Before asking, YES, there are solutions. The desired answers are (x1 = 0.086, y1 = 1.0125) and (x2 = -0.176, y2 = 1.0099), or anything close to that, like (x1 = 0.085, y1 = 1.0126) and (x2 = -0.175, y2 = 1.0097).

How can I find a solution? I appreciate your assistance in advance.

Answer 1

[Sam Chak]

Hi @Faraz j

Fixing the vpasolve() syntax does return some complex-valued solutions, with their real parts close to the values of your user-supplied solutions: (, ) and (, ). Also, check if E1 = 13.7872 and E2 = 11.3728 are rounded to 6 digits of precision. If you use the true values, perhaps you can obtain the real-valued solutions.

% ------------ Parameters --------------------
a       = 0.9;
beta_10 = 3.5668;
k0      = 4.9907;
K2      = 480;

% ------------ 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 --------------------
syms x1 x2 y1 y2

Fn1     = ((cos((beta_10/k0)*y1*v(x1))-cos(y1*v(x1)))/sin(y1*v(x1)))*(sin(pi*x1)/a);
Fn2     = ((cos((beta_10/k0)*y2*v(x2))-cos(y2*v(x2)))/sin(y2*v(x2)))*(sin(pi*x2)/a);

Ya1_G0  = (K2*Fn1^2) / (((K2*Fn1^2)/(g(x1)*h(y1)))+(3.0691-13.7872i));
Ya2_G0  = (K2*Fn2^2) / (((K2*Fn2^2)/(g(x2)*h(y2)))+(-1.9145-11.3728i));

% ------------ Equations --------------------
E1      = imag((K2*Fn1^2)/(g(x1)*h(y1))) == 13.7872; 
E2      = imag((K2*Fn2^2)/(g(x2)*h(y2))) == 11.3728;
E3      = Ya2_G0/Fn2                     == - 2*(Ya1_G0/Fn1);
E4      = Ya1_G0 + Ya2_G0                == 0.5;

E       = [E1,E2,E3,E4];
S       = vpasolve(E, [x1,x2,y1,y2], [a/4 -a/4 0.98 1.02])

S = struct with fields:
    x1: 0.086332800652902918634626942861992 - 0.000075969740964003620819005807052217i
    x2: - 0.17461238106411437018174916235168 + 0.00028056533242226180028405491399577i
    y1: 1.0123659322872014363942075497688 + 0.00069520573785250236695044354710485i
    y2: 1.0094688856718929765239642937614 + 0.0012358492108274508835495455949395i

[Walter Roberson]

If you do not constrain x1 to positive then you get a different solution

% ------------ Parameters --------------------
a = 0.9;
beta_10 = 3.5668;
k0 = 4.9907;
K2 = 480;

% ------------ 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 --------------------
syms x1 x2 y1 y2 real

Fn1 = ((cos((beta_10/k0)*y1*v(x1))-cos(y1*v(x1)))/sin(y1*v(x1)))*(sin(pi*x1)/a);
Fn2 = ((cos((beta_10/k0)*y2*v(x2))-cos(y2*v(x2)))/sin(y2*v(x2)))*(sin(pi*x2)/a);

Ya1_G0 = (K2*Fn1^2) / (((K2*Fn1^2)/(g(x1)*h(y1)))+(3.0691-13.7872i));
Ya2_G0 = (K2*Fn2^2) / (((K2*Fn2^2)/(g(x2)*h(y2)))+(-1.9145-11.3728i));

% ------------ Equations --------------------
E1 = imag((K2*Fn1^2)/(g(x1)*h(y1))) == 13.7872;
E2 = imag((K2*Fn2^2)/(g(x2)*h(y2))) == 11.3728;
E3 = Ya2_G0/Fn2 == -2 * (Ya1_G0/Fn1);
E4 = Ya1_G0 + Ya2_G0 == 0.5;

E = [E1,E2,E3,E4];
S = vpasolve(E,[x1,x2,y1,y2],[0 a/4;-a/4 a/4;0.98 1.02;0.98 1.02])

S = struct with fields:
    x1: 0.086402384311861900570182103043991
    x2: -0.17691470933770734367934875866917
    y1: 1.0125250985379729976135407837632
    y2: 1.0097068286281174578568040509811

S2 = vpasolve(E,[x1,x2,y1,y2],[-a/4 a/4;-a/4 a/4;0.98 1.02;0.98 1.02])

S2 = struct with fields:
    x1: -0.086402384311861900570182103043991
    x2: 0.17691470933770734367934875866917
    y1: 1.0125250985379729976135407837632
    y2: 1.0097068286281174578568040509811

[Sam Chak]

Aha, I forgot to set "assume" the solutions x1, x2, y1, y2 are "real".

[Faraz j]

Thank you so much for your prompt response! Your help is greatly appreciated!

Answer 2

[Torsten]

Try

S = vpasolve(E,[x1,x2,y1,y2],[0.086 -0.176 1.0125 1.0099])

But the solution is complex-valued since your function h is complex-valued.