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您会得到：

作为所有切线的证明：