كيف أُنفّذ طريقة نيوتن رافسون في R ؟
السؤال هو كالتالي :
اريد تعريف متغير (لنفرض x) في خوارزمية نيوتن رافسون لمعادلة انحدار حيث x يمثل متغير مشارك
— Ema. Mas (@EmaMas13) February 27, 2019
طريقة نيوتن رافسون
انا بصراحة ما سمعت عن الطريقة من قبل لذلك بحثت عنها ووجدتها طريقة بسيطة لكنها ذكية. الطريقة تهدف إلى إيحاد قيمة جذر المعادلة بطريقة عددية numerical solution. المعادلة هي كالتالي
\[ x_1 = x_0 - \frac{f(x_0)}{{f}'(x_0)} \]
اولا سوف نقوم بتمثيل هذه المعادلة في R
ثم سنقوم بمتابعت تغير قيمة \(x\) في المعادلة. لكن اولا نحتاج إلى مثال. في صفحة ويكيبيديا هناك مثال لهذه المعادلة
\[
f(x)=x^3-1000x^2+1
\]
سوف نمثلها في الكود كدالة بهذه الطريقة
fx <- function(x){
return(x^3-1000*x^2+1)
}
اشتقاق المعادلة بتطبيق قوانين الإشتقاق سيكون هذه المعادلة \[ {f}'(x)=3x^2-2000x \] سوف نمثلها في الكود كدالة بهذه الطريقة
fx_p <- function(x) {
return(3*x^2-2000*x)
}
عند قراءة طريقة نيوتن رافسون تلاحظ انها حلقة ندور فيها من خلال تعويض المخرجات كمدخلات حتى نصل إلى قيمة الصفر للمعادلة. عندها تكون قيمة المتغير هي جذر المعادلة. هذا النمط من المنطق يعرف بالـ Recursion في علم الحاسب. وهو ان تقوم الدالة بنداء نفسها. لذلك سوف نقوم بتنفيذ الكود على هذا النمط.
newton_raphson <- function(x) {
if (round(fx(x),dp) == 0) return(x)
else {
x_n[(n<<-n+1)] <<- x
x_plus_n <-round((x - (fx(x)/fx_p(x))),20)
return(newton_raphson(x_plus_n))
}
}
في المعادلة ستلاحظ انه يوجد متغيرات جديدة قمت بإضافتها للتحكم بسلوك الدالة فقط وهي
- n
رقم العملية
- x_log
فيكتور يقوم بتسجيل قيم المتغير لكل عملية
- dp
للتحكم بعدد بالمنازل العشرية للنتيجة (إذا كان العدد كبير فسوف تكون النتائج اكثر دقة لكنها قد تسبب في مشكلة stack overflow)
هذه القيم لابد ان تنفذ كل مرة نقوم بعملية الحساب
x_n <- rep(0,300)
n <- 0
dp <- 50
الآن دعنا نرى النتائج
newton_raphson(1)
## [1] 0.03162328
لإستعراض سلوك قيمة \(x\)في المعادلة نقوم بهذه التصوير البياني
library("ggplot2")
df <- data.frame(x_n[which(x_n != 0)],1:10)
names(df) <- c("y", "x")
ggplot(df) +
geom_line(aes(x=x, y = y))+
ylab(expression(x[n]))+
xlab("iterations n")+
theme_bw()
في الختام أرجو ان تكون هذه الخطوات واضحة أن تجيب ولو على جزء كبير من السؤال.
جرب بنفسك
كامل الكود تجده هنا
comments powered by Disqus