найти точку пересечения с помощью scilab - PullRequest
Новая
       25

найти точку пересечения с помощью scilab

1 голос
/ 18 марта 2019

Как найти точки пересечения на графике, показанном ниже, с помощью функции 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)

enter image description here

1 Ответ

0 голосов
/ 19 марта 2019

Ваш оригинальный пост был совершенно нечитабельным и хаотичным. Мне потребовалось время, чтобы отредактировать его и понять, чего вы пытаетесь достичь. Однако я постараюсь вам помочь. Пойдем пошагово:

  1. Я не уверен, почему вы использовали функцию find таким образом. Возможно, вы пытались векторизовать функцию g? Обратите внимание, что Scilab по умолчанию не передает функции через массивы. Вам нужно либо векторизовать их, либо использовать feval для этого. Пожалуйста, прочитайте этот другой ответ Я уже писал ранее. find - векторизованная операция, применяемая к массиву, логическая операция и скаляр, поиск элементов массива, которые удовлетворяют операции. Например, с на страницу find : 
beers = ["Desperados", "Leffe", "Kronenbourg", "Heineken"];
find(beers == "Leffe")

возвращает 2 и 

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

возвращает те элементы массива A, которые меньше 0.4.

  1. Пожалуйста, ознакомьтесь с условиями и, в частности, if, then, elsif, else, end утверждениями. Если вы узнаете это, вы не будете использовать функцию find таким образом. Иногда у вас столько if строк подряд, вместо этого попробуйте использовать 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
  1. Теперь, очевидно, вы хотите найти точки на этой кривой, которые имеют значение 30. Хотя существуют методы для автоматического нахождения этих точек, очень полезно найти правильный диапазон: 
t = [5:50];
plot(t, feval(t, g) - 30)

enter image description here

показывает, что оба решения находятся в диапазоне 20 < x1 < 30 и 40 < x < 50.

  1. Теперь, если мы используем fsolve с правильными начальными значениями, это дает нам хорошие результаты: 
--> deff('[y] = g2(x)', 'y = g(x) - 30');

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

   26.67373
   48.396547
  1. Третий параметр функции fsolve - это Якобин / производная функции g(x). Вы должны либо вычислить производные вышеприведенных полиномов вручную (либо использовать соответствующее символическое программное обеспечение, такое как Maxima), либо определить их как полиномы с помощью функции poly. См. этот урок , например. Затем дифференцируйте их, определяя новую функцию, такую ​​как dgdx.
...