作者:老李
日期:11-1
目标:用解析几何的方法作出椭圆的sino图
需要使用的知识点:
1.符号型变量
2.表达式的求解
3.由矩阵生成灰度图
理论部分:
深刻椭圆表达式
1.X为向量(自变量),X的长度体现为空间的维度。
2.参数 决定了椭圆的中心点位置
3.对 的理解
首先,对 的思考非常有助于对矩阵的理解,是线性代数里非常好的一道题,而且我们还能把这方面的理解拓展到对多维正态分布的理解上。
为一个满秩正定实对称矩阵,它决定了该椭圆的形状。
当
为一个对角矩阵时,它所表现的是一个正规的椭圆(可以理解为单位圆或球往每个轴的方向上的拉伸)
当
不是一个对角矩阵时,它所表现的是一个由正规的椭圆以中心位置旋转之后的椭圆。我们知道,实对称矩阵通过正交变换可以得到一组新的正交基。从这一点的角度上理解,那就是任意一个椭圆(球)在某一个正交基上可以表现为正规椭圆。
我们知道,广义化的正态分布的表达式为:
指数部分里面实际上是一个椭圆的表达式,应为我们知道,协方差矩阵
是正定的,这也完全表达了该概率密度的形状
2.matlab实现:
好吧,我承认我在理论部分扯那么多只是为了掩饰我编程能力的弱小。。。
因为怕麻烦,因为由理论部分我们可以知道,复杂一点的椭圆无非就是对原表达式进行几次线性变换而已,所以我们选用的椭圆为一个非常简单的正规的椭圆。表达式为:
图像为:
我们知到对于直线的表达式可以表示为
我们要计量的是法线角度为0到
之间的情况,而
为直线与原点之间的距离。
在这里插入代码片%定义符号型变量来定义表达式,从而求解
syms x y a b r c;
%求解椭圆与直线联立求解的表达式
[x,y] = solve('x^2/9+y^2/4= 1,a*x+b*y = r');
%将截线长用符号表达出来
c = sqrt((x(1)-x(2))^2+(y(1)^2-y(2))^2);
ezplot('x^2/9+y^2/4 = 1')%作椭圆图
B = zeros(180,length(-3:0.01:3));
v = -3:0.01:3;
for i = 1:180
for k = 1:length(v)
q = (i-1)*pi/180;
d = subs(c,[a;b;r],[sin(q);-cos(q);cos(q)*v(k)]);%以实数替代符号
g = vpa(d);%求值
%滤掉虚数部分,只取有实数解的部分
if abs(imag(g))<eps
B(i,k) = g;
end
end
end
为了避免麻烦,我在写代码的过程中,尽量避免了表达式的推导。
结果,这个程序运行了2000多秒。
然后我优化了一下程序,决定不用subs语句,于是我把表达式先求解了出来。
代码如下:
%(((b*r + a*(4*a^2 + 4*b^2 - r^2)^(1/2))/(a^2 + b^2) - (b*r - a*(4*a^2 + 4*b^2 - r^2)^(1/2))/(a^2 + b^2))^2 + ((r - (b*(b*r + a*(4*a^2 + 4*b^2 - r^2)^(1/2)))/(a^2 + b^2))/a - (r - (b*(b*r - a*(4*a^2 + 4*b^2 - r^2)^(1/2)))/(a^2 + b^2))/a)^2)^(1/2)
B = zeros(180,length(-3:0.01:3));
v = -3:0.01:3;
for i = 1:180
for k = 1:length(v)
q = (i-1)*pi/180+pi/2;
a = sin(q);
b = -cos(q);
r = v(k)*cos(q);
g = (((r - (2*b*(2*b*r - 3*a*(9*a^2 + 4*b^2 - r^2)^(1/2)))/(9*a^2 + 4*b^2))/a - (r - (2*b*(2*b*r + 3*a*(9*a^2 + 4*b^2 - r^2)^(1/2)))/(9*a^2 + 4*b^2))/a)^2 + ((2*(2*b*r - 3*a*(9*a^2 + 4*b^2 - r^2)^(1/2)))/(9*a^2 + 4*b^2) - (2*(2*b*r + 3*a*(9*a^2 + 4*b^2 - r^2)^(1/2)))/(9*a^2 + 4*b^2))^2)^(1/2);
%滤掉虚数部分,只取有实数解的部分
if imag(g)<eps
B(i,k) = g;
end
end
end
%作图
imshow(B,[])
同理,我把表达式变了一下
然后我们得到了新的表达式
`B = zeros(180,length(-3:0.01:3));
v = -3:0.01:3;
for i = 1:180
for k = 1:length(v)
q = (i-1)*pi/180+pi/2;
a = sin(q);
b = -cos(q);
r = v(k)*cos(q);
g = (((r + (b*(a*(a^2 + b^2 - r^2)^(1/2) - b*r))/(a^2 + b^2))/a - (r - (b*(a*(a^2 + b^2 - r^2)^(1/2) + b*r))/(a^2 + b^2))/a)^2 + ((a*(a^2 + b^2 - r^2)^(1/2) + b*r)/(a^2 + b^2) + (a*(a^2 + b^2 - r^2)^(1/2) - b*r)/(a^2 + b^2))^2)^(1/2);
%滤掉虚数部分,只取有实数解的部分
if imag(g)<eps
B(i,k) = g;
end
end
end
%作图
imshow(B,[])
图如下: