Solve a large composite set of equations numerically

I have the following set of equations:

1:

where are known scalar constants, and

2:

3:

where Vis known scalar constant and k is a scalar which I want to find.

4:

5:

I've been informed to find the value of k from the last following two inequalities:

with

and are known constants.

Just to make things little more easier, let

The ultimate problem here is that do not have explicit symbolic equations because I have to compute them numerically since their direct symbolic equations are too complicated and Matlab symbolic engine produces too big output.

So the question is: how to numerically solve inequalities (5:) for the variable k ?

Answer 1

[Andrey Mironenko]

I have finally solved this set symbolically, but this is offtopic, so close it now.

[darova]

Please attach your solution

[Andrey Mironenko]

In the original paper which i've used I have noticed the implicit assumption:

Therefore symbolic evaluation gives:

syms fx fy fz d1 gx gy gz d2 real;
syms x y z real;
syms k positive;
syms s positive;
syms V positive;
syms F G R NC NS;
syms chi_rate_max;
syms alpha positive;

s = sym(1);

f(x,y,z) = fx*x + fy*y + fz*z + d1;
g(x,y,z) = gx*x + gy*y + gz*z + d2;

r(x,y,z) = sqrt(f^2 + g^2);
Nc(x,y,z) = sqrt((f*fx + g*gx)^2 + (f*fy + g*gy)^2 + (f*fz + g*gz)^2);
Ns = sqrt((fy*gz - fz*gy)^2 + (fz*gx - fx*gz)^2 + (fx*gy - fy*gx)^2);

vx_d(x,y,z) = sym(-1)*V*(f*fx + g*gx)*tanh(k*r)/Nc + s*V*(fy*gz - fz*gy)*sech(k*r)/Ns;
vy_d(x,y,z) = sym(-1)*V*(f*fy + g*gy)*tanh(k*r)/Nc + s*V*(fz*gx - fx*gz)*sech(k*r)/Ns;
vz_d(x,y,z) = sym(-1)*V*(f*fz + g*gz)*tanh(k*r)/Nc + s*V*(fx*gy - fy*gx)*sech(k*r)/Ns;

chi_d = atan2(vz_d, vx_d);

B_x = diff(chi_d, x);

B_x = subs(B_x, d1 + fx*x + fy*y + fz*z, F);
B_x = subs(B_x, d2 + gx*x + gy*y + gz*z, G);
B_x = subs(B_x, sqrt(F^2 + G^2), R);
B_x = subs(B_x, R*k, 0);
B_x = subs(B_x, sqrt((F*fx + G*gx)^2 + (F*fy + G*gy)^2 + (F*fz + G*gz)^2), NC);
B_x = subs(B_x, sqrt((fy*gz - fz*gy)^2 + (fz*gx - fx*gz)^2 + (fx*gy - fy*gx)^2), NS);

B_x = simplify(B_x);

B_y = diff(chi_d, y);

B_y = subs(B_y, d1 + fx*x + fy*y + fz*z, F);
B_y = subs(B_y, d2 + gx*x + gy*y + gz*z, G);
B_y = subs(B_y, sqrt(F^2 + G^2), R);
B_y = subs(B_y, R*k, 0);
B_y = subs(B_y, sqrt((F*fx + G*gx)^2 + (F*fy + G*gy)^2 + (F*fz + G*gz)^2), NC);
B_y = subs(B_y, sqrt((fy*gz - fz*gy)^2 + (fz*gx - fx*gz)^2 + (fx*gy - fy*gx)^2), NS);

B_y = simplify(B_y);

B_z = diff(chi_d, z);

B_z = subs(B_z, d1 + fx*x + fy*y + fz*z, F);
B_z = subs(B_z, d2 + gx*x + gy*y + gz*z, G);
B_z = subs(B_z, sqrt(F^2 + G^2), R);
B_z = subs(B_z, R*k, 0);
B_z = subs(B_z, sqrt((F*fx + G*gx)^2 + (F*fy + G*gy)^2 + (F*fz + G*gz)^2), NC);
B_z = subs(B_z, sqrt((fy*gz - fz*gy)^2 + (fz*gx - fx*gz)^2 + (fx*gy - fy*gx)^2), NS);

B_z = simplify(B_z);

ineq1 = sqrt(simplify(B_x^2 + B_z^2)) <= sym(7)*(sym(1)-alpha)*chi_rate_max/(sym(10)*V);
ineq2 = abs(B_y) <= alpha*chi_rate_max/V;

This evaluates into:

ineq1 = (sym)

           _____________________________________________________________________
          /   2                                                                2
         /  NS *((F*fx + G*gx)*(fx*gy - fy*gx) - (F*fz + G*gz)*(fy*gz - fz*gy))
  k*    /   --------------------------------------------------------------------
       /
      /                                   2  2 /               2
    \/                                  NC *R *\(fx*gy - fy*gx)  + (fy*gz - fz*g

  ___________________________________
   /             2                2\
  *\(F*fx + G*gx)  + (F*fz + G*gz) /     chi_rate_max*(7 - 7*alpha)
  ----------------------------------  <= --------------------------
      2                                             10*V
    2\
  y) /
  
  
  ineq2 = (sym)

    |NS*(F*fy + G*gy)*((F*fx + G*gx)*(fx*gy - fy*gx) - (F*fz + G*gz)*(fy*gz - fz
  k*|---------------------------------------------------------------------------
    |                                      NC*R
  ------------------------------------------------------------------------------
                                         2                  2
                          (fx*gy - fy*gx)  + (fy*gz - fz*gy)

  *gy))|
  -----|
       |    alpha*chi_rate_max
  ------ <= ------------------
                    V

From this it is possible to solve inequalities directly with respect to k

[darova]

impossible

[Andrey Mironenko]

Can you please be more specific?

[darova]

Can't believe it's happening

[Andrey Mironenko]

It is if you consider assumption i have pointed above.

Answer 2

[darova]

I think there is no exact solution for k. Try something like this:

% 1
[X,Y,Z] = ndgrid( linspace(-1e5,1e5,20) );
f = fx*X + fy*Y + fz*Z + d1;
% 2
r = sqrt(f.^2+g.^2);
% 3
k = 2;            % assume
vzd = -V*(f*fz+g*gz).*tanh(k*r)./Nc + ...
% 4
gammad = atand2(vzd,vxd);
delta = 2e5/20;
[Bx,By,Bz] = gradient(gammad,delta);
% 5
F1 = sqrt(Bx.^2+Bz.^2) - 7*(1-alpha)*gamma_max/10/V;
F2 = By - alpha*gamma_max/V;
isosurface(X,Y,Z,F1,0)
hold on
isosurface(X,Y,Z,F2,0)
hold off

[Andrey Mironenko]

Just to be clear, are

vzd = 
vxd = 

scalar functions in your answer?

I'm asking because the notation

.* ./

assumes dealing with vectors but

f,g,vzd,vxd

are supposed to be scalar functions.

[darova]

Everything that depend on X,Y,Z are 3D matrices

[Andrey Mironenko]

The code so far:

V = 50;
fx = 0;
fy = 1;
fz = 0;
d1 = -4499;
gx = -1;
gy = 0;
gz = 0;
d2 = -306154;
k = 1e-5;
a = 0.5;
chi_rate_max = 0.4;
% 1
[X,Y,Z] = ndgrid(linspace(-5e6,5e6,20));
f = fx*X + fy*Y + fz*Z + d1;
g = gx*X + gy*Y + gz*Z + d2;
% 2
r = sqrt(f.^2 + g.^2);
Nc = power(((f.*fx + g.*gx).^2 + (f.*fy + g.*gy).^2 + (f.*fz + g.*gz).^2), 0.5);
Ns = power(((fy*gz - fz*gy).^2 + (fz*gx - fx*gz).^2 + (fx*gy - fy*gx).^2), 0.5);
% 3
vzd = -V*(f*fz + g*gz).*tanh(k*r)./Nc + V*(fx*gy - fy*gx).*sech(k*r)./Ns;
vxd = -V*(f*fx + g*gx).*tanh(k*r)./Nc + V*(fy*gz - fz*gy).*sech(k*r)./Ns;
% 4
chid = atan2(vzd,vxd);
delta = 2e5/20;
[Bx,By,Bz] = gradient(chid,delta);
% 5
F1 = sqrt(Bx.^2 + Bz.^2) - 7*(1 - a)*chi_rate_max/10/V;
F2 = By - a*chi_rate_max/V;
isosurface(X,Y,Z,F1,0)
hold on;
isosurface(X,Y,Z,F2,0)
hold off;

gives error:

warning: isosurface: triangulation is empty

warning: called from

isosurface at line 164 column 5

VF at line 31 column 1

error: invalid value for array property "cdata"

error: __go_patch__: unable to create graphics handle

error: called from

__patch__ at line 171 column 7

patch at line 86 column 18

isosurface at line 191 column 10

VF at line 31 column 1