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
采纳的回答
Walter Roberson
2016-5-1
1 个投票
Given those equations, bsfc approximately 1154.643046927194 or approximately 1325.553297436922
Your proposed value 246 does not come at all close to satisfying eqn5
8 个评论
Bibek Tripathi
2016-5-1
Walter Roberson
2016-5-1
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
2016-5-1
Walter Roberson
2016-5-1
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
2016-5-1
Walter Roberson
2016-5-1
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
2016-5-1
Walter Roberson
2016-5-2
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.
更多回答（2 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
