使用scilab找到交叉点

Question

如何使用fsolve函数（来自scilab）在下图中找到交点？

这是我到目前为止所尝试的：

function y=f(x)
    y = 30 + 0 * x;
endfunction


function y= g(x)
    y=zeros(x)
    k1 = find(x >= 5 & x <= 11); 
    if  k1<>[]  then
        y(k1)= -59.535905 +24.763399*x(k1) -3.135727*x(k1)^2+0.1288967*x(k1)^3;
    end;
    k2=find(x >= 11 & x <= 12); 
    if  k2 <> []    then 
        y(k2)=1023.4465 - 270.59543 * x(k2) + 23.715076 * x(k2)^2 - 0.684764 * x(k2)^3; 
    end;
    k3 = find(x >= 12 & x <= 17);    
    if  k3 <> [] then
        y(k3) =-307.31448 + 62.094807 *x(k3) - 4.0091108 * x(k3)^2 + 0.0853523 * x(k3)^3;
    end;
    k4 = find(x >= 17 & x <= 50); 
    if k4 <> [] then 
        y(k4) = 161.42601 - 20.624104 *x(k4) + 0.8567075 * x(k4)^2 - 0.0100559 * x(k4)^3;
    end;
endfunction

t=[5:50];
plot(t, g(t));
plot2d(t, f(t));
deff('res = fct', ['res(1) = f(x)'; 'res(2) = g(x)']);
k1=[5, 45];
xsol1 = fsolve(k1, f, g)

Answer 1

你原来的帖子完全不可读和混乱。我花了一些时间来编辑它并理解你想要实现的目标。但是我会尽力帮助你。让我们一步一步走：

1. 我不知道你为什么用这种方式使用find功能。可能你试图矢量化g功能？请注意，Scilab默认情况下不会在数组上广播函数。你需要矢量化它们或使用feval这样做。请阅读我之前写过的this other answer。 find是一个矢量化操作，应用于数组，布尔运算和标量，找到满足该操作的数组元素。例如来自the find page：

beers = ["Desperados", "Leffe", "Kronenbourg", "Heineken"];
find(beers == "Leffe")

返回2和

A = rand(1, 20);
w = find(A < 0.4)

返回小于A的数组0.4的元素。

2. 请了解条件，特别是if, then, elsif, else, end声明。如果您了解到这一点，就不会以这种方式使用find函数。有时你连续有这么多的ifs，然后尝试用select, case, else, end代替。你的第二个功能可以写成：

function y = g(x)
  if x < 5 | 50 < x then
    error("Out of range");
  elseif x <= 11 then
    y = -59.535905 + 24.763399 * x - 3.135727 * x^2 + 0.1288967 * x^3;
    return;
  elseif x <= 12 then
    y = 1023.4465 - 270.59543 * x + 23.715076 * x^2 - 0.684764 * x^3;
    return;
  elseif x <= 17 then
    y = -307.31448 + 62.094807 * x - 4.0091108 * x^2 + 0.0853523 * x^3;
    return;
  else
    y = 161.42601 - 20.624104 * x + 0.8567075 * x^2 - 0.0100559 * x^3;
  end
endfunction

3. 现在显然你想找到这条曲线上具有30值的点。虽然有方法可以自动绘图找到这些点，但找到合适的范围非常有帮助：

t = [5:50];
plot(t, feval(t, g) - 30)

表明这两种解决方案都在20 < x1 < 30和40 < x < 50的范围内。

4. 现在，如果我们使用具有适当初始值的fsolve，它会给我们带来很好的结果：

--> deff('[y] = g2(x)', 'y = g(x) - 30');

--> fsolve([25; 45], g2)
 ans  =

   26.67373
   48.396547

5. fsolve函数的第三个参数是g(x)函数的Jacobin /导数。您应该手动计算上述多项式的导数（或使用适当的符号软件，如Maxima），或使用poly函数将它们定义为多项式。例如，请参阅this tutorial。然后区分它们，定义一个像dgdx这样的新函数。