solving 5 non linear simultaneous equation

I have a set of 5 non-linear equations. first 4 equations are h,k,a,b --> in terms of variable 'bsfc'. the fifth equation bsfc^2 in terms of h,k,a,b. So i have 5 equations and 5 variables. This programming is basically finding the radius of the rotated ellipse. I have my x and y values too. However I am not getting the right solution.

The right solution is bsfc=246

I have my code here:

syms bsfc h k a b

   y=86;
   x=1717;
   phi=0;
   eqn1= h == 9.514*bsfc;
   eqn2= k == 0.0001*(bsfc+1)^2-0.07*bsfc+19;
   eqn3= a == 10* atan(bsfc/30 - 5.8) -3.5;
   eqn4= b == 0.025 * atan(bsfc/20 - 10.3) -0.022;
   eqn5= bsfc^2 == (((x-h)*cos(phi)+(y-k)*sin(phi))^2/(a^2)+ ((x-h)*sin(phi)-(y-k)*cos(phi))^2/(b^2))

sol=solve([eqn1,eqn2,eqn3,eqn4,eqn5],bsfc,h,k,a,b);
hSol = sol.h
kSol = sol.k
aSol = sol.a
bSol = sol.b
bsfcSol = sol.bsfc

Answer 1

[Walter Roberson]

Given those equations, bsfc approximately 1154.643046927194 or approximately 1325.553297436922

Your proposed value 246 does not come at all close to satisfying eqn5

[Bibek Tripathi]

Thanks for your help, Walter! would changing the following values get me close to bsfc=260?

 y=7.076;
 x=1256.97;
 phi=0.00044;

[Walter Roberson]

With those values, there are solutions at approximately 104.7268201342307, 176.9129085082923, 236.7362969099984, 677.5658592571117

For those x and y, bsfc can be 260 for phi = 0.001964748431506022

[Bibek Tripathi]

That looks interesting, May I please see how you modified the code? Thanks again!

[Walter Roberson]

I left x, y and theta as variables in the equations, and I proceeded backwards from eqn5 towards eqn1, solving for one more variable each time, leaving bsfc for last. After the 5th step, a complicated equation is given for bsfc in terms of x, y, theta.

Then I substituted in the given x and y and theta, and went looking for bsfc that satisfied the equation, which was broadly speaking of the form F(bsfc) = G(bsfc) for functions F and G. Subtract one side from another to get H(bsfc) = F(bsfc) - G(bsfc) and the solutions are where H(bsfc) = 0.

But the equation generates complex values, so I took the imaginary part of H and plotted it over a range of values, looking for the zero crossings and looking for the ranges where it goes 0 for a while. Any solution for bsfc must at least satisfy imag(H(bsfc)) = 0 . Having narrowed down the ranges where solutions are possible in the reals, plot real(H) and imag(H) together on the same graph over nearby ranges and see approximately where real(H) is 0 inside the areas where imag(H) is 0. With those sub-ranges in hand, you can use any root finder on real(H) to find numeric solutions, and then cross-check them against the full H function to be sure that there are no imaginary components.

The determination of phi given bsfc and x and y is much much easier. Go back to the list [eqn1, eqn2, eqn3, eqn4, eqn5] in terms of the variables and x and y and theta, and substitute in the known x and y and the given bsfc . The first 4, [eqn1, eqn2, eqn3, eqn4] then become simple expressions of the form [variable = numeric value]. Substitute all of those into eqn5 to get an equation involving theta. That can be solved for theta. The result is two solutions that are exactly Pi apart, and the whole is periodic on 2*Pi... so effectively the complete set of solutions are spaced Pi apart.

[Bibek Tripathi]

Thanks a lot for your input. Would you be able to share the code somehow? I do not seem to be able to replicate your work. Thanks

[Walter Roberson]

Which direction? Finding bsfc given x and y and theta, or finding theta given x and y and bsfc ?

In the forward direction, finding bsfc, the equation is complex valued, and is impractical or impossible to solve analytically, so it involved taking a bunch of graphs and narrowing down possibilities. There isn't really code for that unless you want to look at a bunch of plot() statements.

The other direction, finding theta, is much easier to express analytically.

[Bibek Tripathi]

I have an array of x and y values. I have to run them for respective bsfc values. So, yes in forward direction. I would be glad to look at your matlab statements and see how you did it. Thanks.

[Walter Roberson]

syms bsfc h k a b x y phi;
Q = @(v) sym(v);
eqn1 = h == Q(9.514)*bsfc;
eqn2 = k == Q(0.1e-3)*(bsfc+1)^2-Q(0.7e-1)*bsfc+19;
eqn3 = a == 10*atan((1/30)*bsfc-Q(5.8))-Q(3.5);
eqn4 = b == Q(0.25e-1)*atan((1/20)*bsfc-Q(10.3))-Q(0.22e-1);
eqn5 = bsfc^2 == ((x-h)*cos(phi)+(y-k)*sin(phi))^2/a^2+((x-h)*sin(phi)-(y-k)*cos(phi))^2/b^2;

sol_hkab = solve([eqn2,eqn3,eqn4,eqn5],[h, k, a, b]);

eqn_bsfc1 = subs(eqn1, [h k a b], [sol_hkab.h(1), sol_hkab.k(1), sol_hkab.a(1), sol_hkab.b(1)]);

eqn_bsfc2 = subs(eqn1, [h k a b], [sol_hkab.h(2), sol_hkab.k(2), sol_hkab.a(2), sol_hkab.b(2)]);

xyt = [Q(1256.97), Q(7.076), Q(0.00044)];

bs = linspace(0,750,2000);

s_bsfc1 = subs(eqn_bsfc1, [x y phi], xyt);
eq1f = matlabFunction(feval(symengine,'lhs',s_bsfc1)-feval(symengine,'rhs',s_bsfc1));
eq1fr = @(x) real(eq1f(x));
eq1fi = @(x) imag(eq1f(x));

s_bsfc2 = subs(eqn_bsfc2, [x y phi], xyt);
eq2f = matlabFunction(feval(symengine,'lhs',s_bsfc2)-feval(symengine,'rhs',s_bsfc2));
eq2fr = @(x) real(eq2f(x));
eq2fi = @(x) imag(eq2f(x));

plot(bs, eq1fi(bs), 'b--', bs, eq2fi(bs), 'g--', bs, eq1fr(bs), 'b-', bs, eq2fr(bs), 'g-');
legend({'imag sol1', 'imag sol2', 'real sol1', 'real sol2'});

and then start zooming in on places where the imaginary is 0, looking for real solutions.

Answer 2

[Alex Sha]

Look this result:

h: -4113.26496133792

bsfc: -432.338129213572

k: 67.8689272162873

a: -18.7135929756326

b: -0.0604868803921316

Answer 3

[Alex Sha]

one more solution:

h: 10985.2739483926

bsfc: 1154.64304691955

k: 71.7260719049847

a: 11.9021369375982

b: 0.0167429176020812