matlab - 找到最大面积(MatLab)

Question

考虑下面的代码(MatLab):

w = 0 : 0.0001 : 9.4978;
a = [1    11    46    95   109    74    24];
b = [-1 3 4 3 1];
mu = 1;
a0 = a(7) ;a1 = a(6) ;a2 = a(5); a3 = a(4) ; a4 = a(3) ; a5 = a(2); a6 = a(1);
b0 = b(5);b1 = b(4);b2 = b(3) ; b3 = b(2); b4 = b(1) ;
De = -a6*w.^6 + a4*w.^4 - a2*w.^2 + a0;
Do = a5*w.^4 - a3*w.^2 + a1;
Ne = b4*w.^4 - b2*w.^2 + b0;
No = -b3*w.^2 + b1;
T = 0.01;
e = real((1i*w).^mu);
f = imag((1i*w).^mu);
A = Ne.*cos(T*w) + w.*No.*sin(T*w);  
B = e.*(Ne.*cos(T*w) + w.*No.*sin(T*w)) - f.*(w.*No.*cos(T*w) - Ne.*sin(T*w));
C = w.*No.*cos(T*w) - Ne.*sin(T*w);
D = e.*(w.*No.*cos(T*w) - Ne.*sin(T*w)) + f.*(Ne.*cos(T*w) + w.*No.*sin(T*w));
Kp = (-De.*D + w.*Do.*B)./(f.*(Ne.^2 + w.^2.*No.^2));
Kd = (-w.*Do.*A + De.*C)./(f.*(Ne.^2 + w.^2.*No.^2));
figure
plot(Kp,Kd)
line([-24 -24],[-2.24 9.813])

通过运行代码，我们得到了这个图:


我想在曲线的指定部分绘制切线( 红色部分 ，w 属于 [0.6342,0.9985] ):


这样做之后，我的目标是找到 最大面积由这条线定义的向内的半平面和由切线产生的所有可能区域之间的曲线(像这样):


在另一点有另一条切线的另一个例子是:



我们可以得出结论，第一个区域比第二个区域大。这种方法应该适用于红色部分的所有点。

我怎样才能通过 MatLab 做到这一点？

我希望我的问题很清楚。任何想法将不胜感激。

Answer 1

这应该以某种方式起作用。

% remove NaNs
Kd(1)=[];
Kp(1)=[];
%%

%exclude non relavant part of original curve
x=Kp;
y=Kd;
exc = 40000;
x(exc:1:end)=[];
y(exc:1:end)=[];

mask = find(x < -9 & x > -19);
xs = x(mask);
ys = y(mask);

L = length(xs)
%%
% determine area of original shape
A_total = polyarea(Kd,Kp);

% pre-allocation    
slope=zeros(L,1)';
inter = slope;
A_part = slope;

for ii = 1:1:L;
    % determine slope for every point
    xslope = xs(ii);
    idx_a = find(xs<xslope,1,'last');
    idx_b = find(xs>xslope,1,'first');
    xa = xs(idx_a);
    xb = xs(idx_b);
    slope(ii) = (ys(idx_b) - ys(idx_a))/(xb - xa);

    % determine slope between current point and any other one
    slopeX = (ys(ii)-y)./(xs(ii)-x);
    % determine intersection points of tangent with rest of curve
    [~,intersection] = min(abs((slopeX)-slope(ii)));
    % index of intersection
    inter(ii)=intersection;
end

% modify curve to get polygon
x_start = x(1);
x_end = x_start;
y_start = y(1);


%finally calculate all single area values A(ii)
for ii = 1:1:L;
    i_inter = inter(ii);
    y_end = y(i_inter) - (x(i_inter)- x_end)*slope(ii);
    x(i_inter+1) = x_end;
    y(i_inter+1) = y_end;
    A_part(ii) = A_total - polyarea( x(1:1:i_inter+1) ,y(1:1:i_inter+1) );
end

当您现在绘图时 A_part在 x你得到:



作为这里所有切线的证明: