Написание функции для поиска корня с использованием метода secant в R - PullRequest
Новая
       48

Написание функции для поиска корня с использованием метода secant в R

0 голосов
/ 17 октября 2018

Метка метода Ньютона-Рафсона неверна, я на самом деле использую метод Секанта;но у меня нет репутации, чтобы создать новый тег.

Я думаю, что я не слишком далеко.У меня есть функция, которую я назвал ISO() для вычисления динамики жидкости.Он нелинейный, поэтому, чтобы найти корень уравнения, мне нужно решить его численно.Я решил использовать метод Secant, так как считаю, что это самый простой метод для того, что я хотел бы сделать.

Я думаю, что у меня возникают проблемы с переназначением первоначальных «догадок», которыеЯ использую, чтобы начать итерацию.Я использую цикл while и подсчитываю количество итераций. 

sec <- function(x){

  Number_Of_Iterations = 0 # Starts the counter

  x1 = x
  x2 = 3*x

  while(x2 - x1 > 0.0001){ # This is how I'm trying to converge the points.

    Number_Of_Iterations = Number_Of_Iterations + 1

    # The secant function to determine a new x:

    x_new = x2 - ISO(x1)*(x2 - x1)/(ISO(x2) - ISO(x1))

    if(x2 - x_new > x_new - x1){ # These were the rules that I set to reassign the inital chosen values.

      x2 = x_new
      x1 = x1

    }else{

      x1 = x_new
      x2 = x2

    }

  }
  m_dot = x_new
  m_dot

}

Когда я запускаю этот код, я получаю число итераций, равное 1, и значение для m_dot равнозначение отличается от того, что достигается при расчете первой итерации вручную.Я проверил, моя функция ISO() возвращает то же значение, что и при вычислении вручную, поэтому базовые функции работают, я просто не думаю, что смогу привести свою функцию sec() в соответствие с написанным мной.

Вся помощь очень ценится.

Ответы [ 2 ]

0 голосов
/ 18 октября 2018

К сожалению, моя проблема была не просто проблемой кодирования;это было проблемой при использовании метода Secant.Мои начальные оценки для начала итераций заставили функцию анализировать неверный диапазон.

Совет, данный LutzL, очень долго отвечал на мой вопрос, abs() было необходимо, и его (или ее) предложение убрать функции из секущей формулы, вероятно, будет использовать гораздо меньше вычисленийpower.

Вот мой рабочий код: 

sec <- function(x){

  x1 = x
  x2 = x/10 # I have changed the initial values so that x1 is the larger of the two initial estimates.

  while(abs(x2 - x1) > 0.0000001){ # This is how I'm trying to converge the points.

    # The secant function to determine a new x:

    x_new = x2 - square(x2)*(x2 - x1)/(square(x2) - square(x1))

    # I have changed from ISO(), to a simple quadratic formula, I named square(). This secant function works for all of
    # the functions I have since tested it on.

    if( abs(x2 - x_new) > abs(x_new - x1)){

      x2 = x_new

    }else{

      x1 = x_new

    }

  }

  m_dot = x_new
  m_dot

}

Причина, по которой я решил написать свою собственную секущую формулу вместо использования встроенной функции из пакета, заключалась в том, что эту функцию очень просто создатьи я думаю, что более полезно использовать мой собственный код, кроме как полагаться на пакеты, написанные другими.

Я надеюсь, что это поможет другим, новичкам в кодировании, и большое спасибо тем, кто предоставил ответы, особенно LutzLкак и без их предложений, я не думаю, что смог бы решить эту проблему так же быстро, как я.

0 голосов
/ 18 октября 2018

Я могу догадаться, что произошло: значения функций в начальных точках имеют одинаковый знак, так что секущий корень x_new находится вне интервала [x1, x2].Затем с помощью логики повторного назначения вы получите это на следующей итерации x2<x1, так что при проверке длины интервала разница x2-x1 будет отрицательной и, следовательно, будет меньше положительного порога.Таким образом, выход после одной итерации без поиска корня.

Наиболее вероятное намерение состояло в том, что все сравнения были для длины, то есть абсолютной величины различий.Таким образом, измените на 

 while( abs(x2 - x1) > 0.0001){ # This is how I'm trying to converge the points.
    ...
    if( abs(x2 - x_new) > abs(x_new - x1) ){ # These were the rules that I set to reassign the inital chosen values.

Кстати: по вашему описанию, функция ISO кажется дорогостоящей.Сохраните значения их оценки один раз и используйте их снова, где это необходимо, в основном в 

x_new = x2 - y2*(x2 - x1)/(y2 - y1)

(обратите внимание, что я исправил опечатку с y1 до y2, вы решаете равенство наклонов (x2-x_new)/(y2-0) = (x2 - x1)/(y2 - y1), чтобы найтисекущий корень.)

Таким образом, установив y1=ISO(x1); y2=ISO(x2); в начале и внутри цикла, добавьте к переназначению значений x вычисление значения y. 

   x2 = x_new; y2 = ISO(x_new);

соответственно 

   x1 = x_new; y1 = ISO(x_new);

Таким образом, для любого из x используемых значений функция ISO оценивается ровно один раз.

...