[已答复] matlab绘制带有限制条件的隐函数

[复制链接]
1940068 发表于 2021-9-3 09:06:51
10 财富积分

本帖最后由 1940068 于 2021-9-8 15:42 编辑

为什么红色线(代码在indx1)的下部分变成斜线,而不是和黑线(代码在indx2)一起掩盖青线。(-280)是函数,indx1内的是条件,有什么方法可以让黑色线和红色线一起掩盖青色线?
f1=ezplot(' -280.3.*y-1.29.*x-426.05.*0.001+log10(y)+log10(0.001)+350137.*0.001.^2-20050.*y.^2-0.01.*x.^2-150.*x.*y+401462.*y.*0.001+72.05.*x.*0.001+6.79',[10^-4,10^1,10^-5,10^-1]);
set(gca,'xscale','log','yscale','log');

hold on
f1=ezplot(' -280.3.*y-1.29.*x-426.05.*0.001+log10(y)+log10(0.001)+350137.*0.001.^2-20050.*y.^2-0.01.*x.^2-150.*x.*y+401462.*y.*0.001+72.05.*x.*0.001+6.79',[10^-4,10^1,10^-5,10^-1]);
set(gca,'xscale','log','yscale','log');
xy=get(f1,'ContourMatrix');
x1=xy(1,2:end);
y1=xy(2,2:end);
indx1=(find(-1124.46.*y1-4.72.*x1-430.52.*0.001-80200.*y1.^2-0.04.*x1.^2+350216.64.*0.001.^2+178826.*y1.*0.001+978.14.*x1.*0.001-600.*x1.*y1+log10(y1)+2.*log10(x1)+4.*log10(0.001)+20.54<=0));
indx2=(find(-1124.46.*y1-4.72.*x1-430.52.*0.001-80200.*y1.^2-0.04.*x1.^2+350216.64.*0.001.^2+178826.*y1.*0.001+978.14.*x1.*0.001-600.*x1.*y1+log10(y1)+2.*log10(x1)+4.*log10(0.001)+20.54>=0));
plot(x1(indx1),y1(indx1),'r',x1(indx2),y1(indx2),'k');
untitled.jpg

2 条回复


TouAkira 发表于 2021-9-3 15:07:33
自己多琢磨琢磨啊,你把几个曲面、曲线全画出来再看,不就明白了嘛
  1. x = 10 .^ [ -4 : 0.25 : 1 ];
  2. y = 10 .^ [ -4 : 0.1 : -1 ];
  3. [ X, Y ] = meshgrid( x, y );
  4. Z = @( x, y ) -280.3.*y-1.29.*x-426.05.*0.001+log10(y)+log10(0.001)+350137.*0.001.^2-20050.*y.^2-0.01.*x.^2-150.*x.*y+401462.*y.*0.001+72.05.*x.*0.001+6.79;
  5. Z_1 = Z( X, Y );
  6. ConstraintFunc = @( x1, y1 ) -1124.46.*y1-4.72.*x1-430.52.*0.001-80200.*y1.^2-0.04.*x1.^2+350216.64.*0.001.^2+178826.*y1.*0.001+978.14.*x1.*0.001-600.*x1.*y1+log10(y1)+2.*log10(x1)+4.*log10(0.001)+20.54;
  7. Z_2 = ConstraintFunc( X, Y );
  8. f1=ezplot(' -280.3.*y-1.29.*x-426.05.*0.001+log10(y)+log10(0.001)+350137.*0.001.^2-20050.*y.^2-0.01.*x.^2-150.*x.*y+401462.*y.*0.001+72.05.*x.*0.001+6.79',[10^-4,10^1,10^-5,10^-1]);
  9. set(gca,'xscale','log','yscale','log');
  10. xy=get(f1,'ContourMatrix');
  11. x1=xy(1,2:end);
  12. y1=xy(2,2:end); close all;
  13. indx1=(find(-1124.46.*y1-4.72.*x1-430.52.*0.001-80200.*y1.^2-0.04.*x1.^2+350216.64.*0.001.^2+178826.*y1.*0.001+978.14.*x1.*0.001-600.*x1.*y1+log10(y1)+2.*log10(x1)+4.*log10(0.001)+20.54<0));
  14. indx2=(find(-1124.46.*y1-4.72.*x1-430.52.*0.001-80200.*y1.^2-0.04.*x1.^2+350216.64.*0.001.^2+178826.*y1.*0.001+978.14.*x1.*0.001-600.*x1.*y1+log10(y1)+2.*log10(x1)+4.*log10(0.001)+20.54>=0));

  15. s1 = surf( X, Y, Z_1, 'FaceAlpha', 0.2, 'LineStyle', '-', 'LineWidth', 0.05 );
  16. hold on; view( 2 ); grid on; zlim( [ -5, 5 ] );
  17. s2 = surf( X, Y, Z_2, 'FaceColor', 'b', 'FaceAlpha', 0.2, 'LineStyle', '-', 'LineWidth', 0.05 );
  18. set(  gca, 'xscale', 'log', 'yscale', 'log' );
  19. c1 = contour( X, Y, Z_1, [ 0, 0 ], 'm-', 'LineWidth', 2 );
  20. c2 = contour( X, Y, Z_2, [ 0, 0 ], 'b-', 'LineWidth', 2 );   
  21. plot3( x1( indx1 ), y1( indx1 ), zeros( size( x1(indx1) ) ), 'r:', ...
  22.     x1( indx2 ), y1( indx2 ), zeros( size( x1(indx2) ) ), 'k-.', 'LineWidth', 3 );
复制代码

s1是你用ezplot绘制的函数的曲面,f1实际上应当是该曲面高度为零的等高线,即c1,但由于插值精度之类问题,f1与c1之间可能会有些许偏差。
s2是你的判断条件对应曲面,按与零的关系区分indx1与indx2,所以c2绘制了它的高度为零的等高线。
然后就很明显了,c1 c2是交叉的。indx2对应的部分位于c2内部,indx1对应曲线位于其外部。
c1曲线图形如下
indx1, c1, indx1, indx1, ...          indx1
                                                c2
                                              indx2
indx1, c2, indx2, indx2, ...          indx2
注意它的左下部分,是由 indx1, c2, indx2 这样的插值点相连的,如果插值不够密集的话,比如
indx1对应横坐标为0, c2对应横坐标为50, indx2对应横坐标为100
那么当你非要区分indx1与indx2时,右上角的indx1直接与左下角的indx1相连接,c2内部的indx2自己连接,而横坐标0到100之间,本来在c1曲线上经插值点相连的曲线,这时候就什么都没有了,看起来当然就“断”了。

1940068 发表于 2021-9-3 15:27:42
TouAkira 发表于 2021-9-3 15:07
自己多琢磨琢磨啊,你把几个曲面、曲线全画出来再看,不就明白了嘛

s1是你用ezplot绘制的函数的曲面,f1实 ...

那有什么办法可以处理吗?我主要是想让红色的线下部分掩盖那段青线

1940068 发表于 2021-9-5 11:16:56
我主要是想画镁铝尖晶石、MgO和氧化铝的优势区域图,遇到此步就卡住了,求大佬帮帮忙
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表