椭圆的交点求解

ivylife

看到这个论坛高手如云，想请教一下 如何解一个二元二次方程组。其实是两个椭圆的交点，要求求出4个解。想得到的是其中的复数解。我用solve解的，带到方程里发现根本不等于0，用fsolve只能得到1个解。

[ 本帖最后由 ChaChing 于 2009-2-23 16:08 编辑 ]

ChaChing

两个椭圆的交点会是复数解!?
建议LZ直接给出完整程序, 方便有兴趣者试!

xukai871105

使用fsolve命令的话只能有一个解

你要想得到所有解的话,可以适当改变初始值

ivylife

具体问题是这样的。（1 ，i, 0）(1,-i,0)是一对圆环点，任何一个圆都经过这两个点。当两个圆经历了投影变换，就变成了两个椭圆，我需要把这两个椭圆的交点解出来，以求出圆环点在这个投影变换中的像。理论上来说是四个解，其中一对是圆环点的像。二元二次方程组的系数有些复杂， X=(x,y,1)
e1= [   0.0192   -0.0016   -0.1880;   -0.0016    0.0255    0.0012;   -0.1880    0.0012    1.8330];
e2= [   0.0188   -0.0014   -0.1742;   -0.0014    0.0256   -0.0081;   -0.1742   -0.0081    1.6216];
这是一个系数矩阵。
方程是 X.'*E1*X=0;   X.'*E2*X=0;
系数矩阵和多项式的转化关系是
% a11x^2+a22y^2+2a13x+2a23y+a33=0

% (a11 a12 a13)
% (a12 a22 a23)  =A
% (a13 a23 a33)
定义的函数是。
function F = myfun(x,a1,a2,a3,a4,a5,b1,b2,b3,b4,b5)
F = [a1*x(1)^2+a2*x(2)^2+a3*x(1)+a4*x(2)+a5;
      b1*x(1)^2+b2*x(2)^2+b3*x(1)+b4*x(2)+b5];
具体算的时候是：
a1=e1(1,1); a2=e1(2,2); a3=e1(1,3)*2; a4=e1(2,3)*2; a5=e1(3,3);
b1=e2(1,1); b2=e2(2,2); b3=e2(1,3)*2; b4=e2(2,3)*2; b5=e2(3,3);
x0=[0.5 0.5]; options=optimset('Display','iter');   % Option to display output
[x,fval] = fsolve(@myfun,x0,options,a1,a2,a3,a4,a5,b1,b2,b3,b4,b5)
disp(x);

[ 本帖最后由 ChaChing 于 2009-2-23 20:49 编辑 ]

ivylife

我已经在这个问题上纠缠了三个星期了，高手们，帮帮我吧。

ChaChing

可能专业不同及个人水平有限! 看了半天, 还是不清楚楼主的问题!
圆环点(1,i,0)(1,-i,0)中, i是虚数吗? 不懂! 不是空间上的点?
投影变换是否纯粹投影至X-Y平面? 为何一定有四个解?
不知是否误会了楼主的问题!

ivylife

i 是虚数单位。是在射影空间里的点。其实就是相机成像过程里的问题。

投影变换是从三维世界空间，变换到二维图像空间。
其实这些都不重要。关键是如何能解出这个方程组来。一个二元二次方程组，也就是两个椭圆的曲线，解方程，为什么就这么难？？？？

ChaChing

对不起, 不懂就喜欢问! 其实喜欢问清楚尚有其他目的, 工程问题不见得仅有一种路径
fsolve一次仅能解一个解, 若许LZ知道答案在那儿附近, 试试不同的初值x0看看!
不然仅可试试solve解了!
不懂"带到方程里发现根本不等于0"的涵义? 差很多吗?

ivylife

晕，说什么对不起呢。我是对这个东西很苦恼，欢迎大家问我啊。呵呵。
我也很奇怪，差别真的很大。solve能解出四个解来，带入方程，结果不是0.

ivylife

原来是我把椭圆的xy相给漏了。。。。。太粗心大意了。。。
用solve，解决了。。都解出来了。。。
非常感谢所有回复的人。尤其感谢chaching的提示。好运。。。

[ 本帖最后由 ivylife 于 2009-2-23 12:42 编辑 ]

ChaChing

汗! 根本没提示什麽!

ChaChing

LS有空时整理下, 发个最後程序来! 作个ending!
让後来者可以学习! 考虑下!