برنامه نویسی، ریاضی ۱۲۷۰۹ بازدید

در محاسبات عددی یا همان آنالیز عددی (Numerical Analysis)، «روش‌های رانگ کوتا» (Runge–Kutta Methods) (روش‌های رونگه کوتا نیز به آن گفته می‌شود) خانواده‌ای از روش‌های تکرار شدنی «صریح و ضمنی» (Explicit and Implicit Methods) هستند. یکی از این روش‌ها، یعنی رانگ کوتای مرتبه اول با عنوان «روش اویلر» نیز شناخته شده و از محبوبیت قابل توجهی برخوردار است. روش‌های رانگ کوتا برای «گسسته‌سازی زمانی» (Temporal Discretization) جهت تخمین راه حل‌های «معادلات دیفرانسیل معمولی» (Ordinary Differential Equations) مورد استفاده قرار می‌گیرند. این روش‌ها در حدود سال ۱۹۰۰ توسط دو ریاضی‌دان آلمانی به نام‌های «کارل رانگ» (Carl Runge) و «ویلهلم کوتا» (Wilhelm Kutta) ساخته شده است. شناخته شده‌ترین روش از خانواده روش‌های رانگ کوتا، RK رانگ کوتای کلاسیک و یا روش رانگ کوتای مرتبه چهارم است که به طور خلاصه به آن روش رانگ کوتا گفته می‌شود. در این مطلب، به روش رانگ کوتا پرداخته شده و کد پیاده‌سازی آن در زبان‌های گوناگون ارائه شده است. همچنین، برای درک بهتر مطلب، مثال هایی از حل معادلات دیفرانسیل با روش رانگ کوتا ارائه شده و پیاده‌سازی آن در زبان برنامه‌نویسی «متلب» (MATLAB) انجام شده است.

محتوای این مطلب جهت یادگیری بهتر و سریع‌تر آن، در انتهای متن به صورت ویدیویی نیز ارائه شده است.

برای مشاهده ویدیوها کلیک کنید.

روش رانگ کوتا مرتبه چهارم

ورودی‌های زیر وجود دارند:

  • یک معادلات دیفرانسیل معمولی که مقدار dy/dx را به شکل x و y تعریف می‌کند.
  • مقدار اولیه y (یعنی (0)y)

بدین ترتیب، معادله زیر وجود دارد.

اکنون، باید مقدار تابع ناشناخته y در نقطه x را به دست آورد. روش رانگ کوتا مقدار تقریبی y برای یک x داده شده را پیدا می‌کند. تنها، معادلات دیفرانسیل معمولی مرتبه اول را می‌توان با استفاده از روش رانگ کوتا مرتبه چهارم حل کرد. در زیر، فرمول مورد استفاده برای محاسبه مقدار بعدی yn+1 از مقدار پیشین yn به دست می‌آید. مقدار n برابر $$0, 1, 2, 3, ….(x – x_0)/h$$ است. در اینجا، h طول هر گام و xn+1 = x0 + h است. طول گام کمتر به معنای صحت بیشتر است.

فرمول بالا اساسا مقدار بعدی yn+1 را با استفاده از yn کنونی، همراه با میانگین چهار افزایش محاسبه می‌کند.

  • k1 افزایش بر پایه شیب در آغاز بازه با استفاده از y است.
  • k2 افزایش بر پایه شیب در نقطه میانی بازه با استفاده از y + hk1/2 است.
  • k3 افزایش بر پایه نقطه میانی با استفاده از y + hk2/2 است.
  • k4 افزایش بر پایه شیب در پایان بازه با استفاده از y + hk3 است.

متد بیان شده در بالا، روش مرتبه چهارم است. بدین معنا که خطای برش محلی از مرتبه (O(h5 است؛ در حالی که خطای کلی تجمیع شده از مرتبه (O(h4 است. در ادامه، پیاده‌سازی این متد در زبان‌های برنامه‌نویسی «پایتون» (Python)، «جاوا» (Java)، «سی» (C) و «سی‌پلاس‌پلاس» (++C)، «پی‌اچ‌پی» (PHP) و «سی‌شارپ» (#C) ارائه شده است.

کد روش رانگ کوتا در پایتون

# Python program to implement Runge Kutta method 
# A sample differential equation "dy / dx = (x - y)/2" 
def dydx(x, y): 
    return ((x - y)/2) 
  
# Finds value of y for a given x using step size h 
# and initial value y0 at x0. 
def rungeKutta(x0, y0, x, h): 
    # Count number of iterations using step size or 
    # step height h 
    n = (int)((x - x0)/h)  
    # Iterate for number of iterations 
    y = y0 
    for i in range(1, n + 1): 
        "Apply Runge Kutta Formulas to find next value of y"
        k1 = h * dydx(x0, y) 
        k2 = h * dydx(x0 + 0.5 * h, y + 0.5 * k1) 
        k3 = h * dydx(x0 + 0.5 * h, y + 0.5 * k2) 
        k4 = h * dydx(x0 + h, y + k3) 
  
        # Update next value of y 
        y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4) 
  
        # Update next value of x 
        x0 = x0 + h 
    return y 
  
# Driver method 
x0 = 0
y = 1
x = 2
h = 0.2
print 'The value of y at x is:', rungeKutta(x0, y, x, h) 
  
# This code is contributed by Prateek Bhindwar

کد روش رانگ کوتا در جاوا

// Java program to implement Runge Kutta method 
import java.io.*; 
class differential 
{ 
    double dydx(double x, double y) 
    { 
        return ((x - y) / 2); 
    } 
      
    // Finds value of y for a given x using step size h 
    // and initial value y0 at x0. 
    double rungeKutta(double x0, double y0, double x, double h) 
    { 
        differential d1 = new differential(); 
        // Count number of iterations using step size or 
        // step height h 
        int n = (int)((x - x0) / h); 
  
        double k1, k2, k3, k4, k5; 
  
        // Iterate for number of iterations 
        double y = y0; 
        for (int i = 1; i <= n; i++)  
        { 
            // Apply Runge Kutta Formulas to find 
            // next value of y 
            k1 = h * (d1.dydx(x0, y)); 
            k2 = h * (d1.dydx(x0 + 0.5 * h, y + 0.5 * k1)); 
            k3 = h * (d1.dydx(x0 + 0.5 * h, y + 0.5 * k2)); 
            k4 = h * (d1.dydx(x0 + h, y + k3)); 
  
            // Update next value of y 
            y = y + (1.0 / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4); 
              
            // Update next value of x 
            x0 = x0 + h; 
        } 
        return y; 
    } 
      
    public static void main(String args[]) 
    { 
        differential d2 = new differential(); 
        double x0 = 0, y = 1, x = 2, h = 0.2; 
          
        System.out.println("\nThe value of y at x is : " 
                    + d2.rungeKutta(x0, y, x, h)); 
    } 
} 
  
// This code is contributed by Prateek Bhindwar

کد روش رانگ کوتا در C++/C

// C program to implement Runge Kutta method 
#include<stdio.h> 
  
// A sample differential equation "dy/dx = (x - y)/2" 
float dydx(float x, float y) 
{ 
    return((x - y)/2); 
} 
  
// Finds value of y for a given x using step size h 
// and initial value y0 at x0. 
float rungeKutta(float x0, float y0, float x, float h) 
{ 
    // Count number of iterations using step size or 
    // step height h 
    int n = (int)((x - x0) / h); 
  
    float k1, k2, k3, k4, k5; 
  
    // Iterate for number of iterations 
    float y = y0; 
    for (int i=1; i<=n; i++) 
    { 
        // Apply Runge Kutta Formulas to find 
        // next value of y 
        k1 = h*dydx(x0, y); 
        k2 = h*dydx(x0 + 0.5*h, y + 0.5*k1); 
        k3 = h*dydx(x0 + 0.5*h, y + 0.5*k2); 
        k4 = h*dydx(x0 + h, y + k3); 
  
        // Update next value of y 
        y = y + (1.0/6.0)*(k1 + 2*k2 + 2*k3 + k4);; 
  
        // Update next value of x 
        x0 = x0 + h; 
    } 
  
    return y; 
} 
  
// Driver method 
int main() 
{ 
    float x0 = 0, y = 1, x = 2, h = 0.2; 
    printf("\nThe value of y at x is : %f", 
            rungeKutta(x0, y, x, h)); 
    return 0; 
}

کد روش رانگ کوتا در PHP

<?php 
// PHP program to implement  
// Runge Kutta method 
  
// A sample differential equation 
// "dy/dx = (x - y)/2" 
function dydx($x, $y) 
{ 
    return(($x - $y) / 2); 
} 
  
// Finds value of y for a 
// given x using step size h 
// and initial value y0 at x0. 
function rungeKutta($x0, $y0, $x, $h) 
{ 
      
    // Count number of iterations 
    // using step size or step  
    // height h 
    $n = (($x - $x0) / $h); 
  
    $k1; $k2; $k3; $k4; $k5; 
  
    // Iterate for number  
    // of iterations 
    $y = $y0; 
    for($i = 1; $i <= $n; $i++) 
    { 
          
        // Apply Runge Kutta 
        // Formulas to find 
        // next value of y 
        $k1 = $h * dydx($x0, $y); 
        $k2 = $h * dydx($x0 + 0.5 * $h, 
                        $y + 0.5 * $k1); 
        $k3 = $h * dydx($x0 + 0.5 * $h,  
                        $y + 0.5 * $k2); 
        $k4 = $h * dydx($x0 + $h, $y + $k3); 
  
        // Update next value of y 
        $y = $y + (1.0 / 6.0) * ($k1 + 2 *  
                    $k2 + 2 * $k3 + $k4);; 
  
        // Update next value of x 
        $x0 = $x0 + $h; 
    } 
  
    return $y; 
} 
  
    // Driver method 
    $x0 = 0; 
    $y = 1;  
    $x = 2;  
    $h = 0.2; 
    echo "The value of y at x is : ", 
          rungeKutta($x0, $y, $x, $h); 
  
// This code is contributed by anuj_67. 
?>

کد روش رانگ کوتا در #C

// C# program to implement Runge 
// Kutta method 
using System; 
  
class GFG { 
      
    static double dydx(double x, double y) 
    { 
        return ((x - y) / 2); 
    } 
      
    // Finds value of y for a given x 
    // using step size h and initial 
    // value y0 at x0. 
    static double rungeKutta(double x0, 
                double y0, double x, double h) 
    { 
      
        // Count number of iterations using 
        // step size or step height h 
        int n = (int)((x - x0) / h); 
  
        double k1, k2, k3, k4; 
  
        // Iterate for number of iterations 
        double y = y0; 
          
        for (int i = 1; i <= n; i++)  
        { 
              
            // Apply Runge Kutta Formulas 
            // to find next value of y 
            k1 = h * (dydx(x0, y)); 
              
            k2 = h * (dydx(x0 + 0.5 * h, 
                             y + 0.5 * k1)); 
                               
            k3 = h * (dydx(x0 + 0.5 * h, 
                            y + 0.5 * k2)); 
                              
            k4 = h * (dydx(x0 + h, y + k3)); 
  
            // Update next value of y 
            y = y + (1.0 / 6.0) * (k1 + 2  
                       * k2 + 2 * k3 + k4); 
              
            // Update next value of x 
            x0 = x0 + h; 
        } 
          
        return y; 
    } 
      
    // Driver code 
    public static void Main() 
    { 
          
        double x0 = 0, y = 1, x = 2, h = 0.2; 
          
        Console.WriteLine("\nThe value of y"
                             + " at x is : "
                 + rungeKutta(x0, y, x, h)); 
    } 
} 
  
// This code is contributed by Sam007.

خروجی

در کد بالا، به ازای y=1، x=2، x0=0 و h=0.2، خروجی به صورت زیر خواهد بود. باید به این نکته نیز توجه کرد که پیچیدگی روش ارائه شده در بالا برابر با (O(n است که در آن، n برابر با x-x0)/h) است.

The value of y at x is : 1.103639

حل معادله دیفرانسیل با روش رانگ کوتا مرتبه چهارم و پیاده‌سازی در متلب

در ادامه، دو مثال ارائه شده است که طی آن‌ها، یک معادله دیفرانسیل با روش رانگ کوتا مرتبه چهارم حل و پیاده‌سازی آن در زبان برنامه‌‌نویسی متلب انجام شده است.

مثال اول

برنامه‌ای بنویسید که معادله دیفرانسیل زیر را با طول گام دلخواه حل کند.

روش رانگ کوتا مرتبه چهارم (RK4) -- از صفر تا صد

حل معادله دیفرانسیل بالا با روش رانگ کوتا مرتبه چهارم و پیاده‌سازی آن با زبان برنامه‌نویسی متلب انجام شده است.

clc;
clear;
close all;
%%
h = 0.1 ; 
t = 0:h:2 ; 
x = t ; 
Nt = numel(t) ; 
y = zeros(1 , Nt) ; y(1) = 0.5 ; 
f = @(t , y) y - t^2 +1 ;

for i = 1:Nt-1
    K1 = h * f(x(i) , y(i)) ; 
    K2 = h * f(x(i)+h/2 , y(i)+K1/2) ; 
    K3 = h * f(x(i)+h/2 , y(i)+K2/2) ; 
    K4 = h * f(x(i)+h , y(i)+K3) ; 
    y(i+1) = y(i) + 1/6*(K1 + 2*K2 + 2*K3+K4) ;
end

yr = (1+t).^2-0.5.*exp(t); 


figure(1) ; 
subplot(211);
plot(t , yr , t , y , 'LineWidth' , 3) ;
xlabel('t'); ylabel('y') ; 
grid on
legend('Real Y' , 'App Y');

subplot(212);
plot(t , yr-y , 'LineWidth' , 3) ;
xlabel('t'); ylabel('error') ; 
grid on
legend('Error');

مثال دوم

برنامه‌ای بنویسید که معادله دیفرانسیل زیر را با طول گام دلخواه حل کند.

روش رانگ کوتا مرتبه چهارم (RK4) -- از صفر تا صد

حل معادله دیفرانسیل بالا با روش رانگ کوتا مرتبه چهارم و پیاده‌سازی آن با زبان برنامه‌نویسی متلب انجام شده است.

clc;
clear;
close all;
%%
h = 0.1; 
x = 0:h:5 ; 
Nt = numel(x) ; 
y = zeros(2 , Nt) ; y(: , 1) = [0.5;0] ; 
f = @(x , y) [y(2) ; x*y(1)];

for i = 1:Nt-1
    K1 = h * f(x(i) , y(: , i)) ; 
    K2 = h * f(x(i)+h/2 , y(: , i)+K1/2) ; 
    K3 = h * f(x(i)+h/2 , y(: , i)+K2/2) ; 
    K4 = h * f(x(i)+h , y(: , i)+K3) ; 
    y(: , i+1) = y(:  , i) + 1/6*(K1 + 2*K2 + 2*K3+K4) ;
end

figure(1) ; 
plot(x , y , 'LineWidth' , 3) ;
xlabel('x'); ylabel('y') ; 
grid on
legend('Y1 = z(x)' , 'Y2 = dot Z(x)');

اگر نوشته بالا برای شما مفید بوده است، آموزش‌های زیر نیز به شما پیشنهاد می‌شوند:

^^

فیلم‌ های آموزش روش رانگ کوتا مرتبه چهارم (RK4) — از صفر تا صد (+ دانلود فیلم آموزش گام به گام)

فیلم آموزشی روش رانگ‌ کوتای مرتبه چهار

دانلود ویدیو

فیلم آموزشی حل مثال از روش رانگ‌‌کوتای مرتبه چهار

دانلود ویدیو

فیلم آموزشی کدنویسی روش رانگ‌ کوتای مرتبه چهار در متلب

دانلود ویدیو
بر اساس رای ۳۸ نفر
آیا این مطلب برای شما مفید بود؟
شما قبلا رای داده‌اید!
اگر بازخوردی درباره این مطلب دارید یا پرسشی دارید که بدون پاسخ مانده است، آن را از طریق بخش نظرات مطرح کنید.

«الهام حصارکی»، فارغ‌التحصیل مقطع کارشناسی ارشد مهندسی فناوری اطلاعات، گرایش سیستم‌های اطلاعات مدیریت است. او در زمینه هوش مصنوعی و داده‌کاوی، به ویژه تحلیل شبکه‌های اجتماعی، فعالیت می‌کند.

12 نظر در “روش رانگ کوتا مرتبه چهارم (RK4) — از صفر تا صد (+ دانلود فیلم آموزش گام به گام)

  • با سلام. من یک سوالی داشتم. برای معادله دیفرانسیل مرتبه دو باید چکار کنم؟ یعنی معادله ای به شکل زیر:
    d^2(x)/dt^2+(c/m)*dx/dt+(k/m)*x
    این یک معادله ارتعاشی می باشد و برابر با صفر می باشد. ممنون میشم اگر راهنمای ی بفرمایید.

  • با سلام و خسته نباشید.
    ببخشید توی مثال اول اگر (y(t بشه (۳۰-y(t یا هر عدد به جای ۳۰ اونوقت باید چکار کنیم.

  • سلام کاش برا معادله دیفرانسیل مرتبه سوم هم میزاشتین خیلی بدردم میخوره
    و برای شرایط مرزی اگه ۰=(بینهایت)”y بجای ۰=(۰)”y بدن چی میشه؟
    ممنون میشم بگین

    1. سلام. برای یادگیری روش حل معادله دیفرانسیل مرتبه سوم، پیشنهاد می‌کنیم دو مطلب «روش اویلر — از صفر تا صد (+ دانلود فیلم آموزش گام به گام)» و «روش رانگ کوتای مرتبه دوم — از صفر تا صد» را مطالعه کنید.
      از همراهی‌تان با مجله فرادرس خوشحالیم.

  • سعید خانی — says: ۳ مهر، ۱۳۹۸ در ۴:۲۴ ب٫ظ

    سلام.
    کد متلب نمونه مثال شما رو اماده کردم .
    انشالله که مفدی واقع بشه.

    %% Runge Kutta method

    clc
    clear all
    close all
    %% Driver method
    x0 = 0;
    y = 1;
    x = 2;
    h = 0.2;
    y = rk4(x0, y, x, h);
    xi = linspace(x0,x,(x-x0)/h);
    display(sprintf(‘The value of y at x is : %f ‘, y(end)));

    plot(xi,y);
    xlabel(‘x’)
    ylabel(‘f(x,y)’)
    %% A sample differential equation
    %% “dy/dx = (x – y)/2”
    function f = dydx(x,y)
    f = (x-y) / 2;
    end

    %% Finds value of y for a
    %% given x using step size h
    %% and initial value y0 at x0.
    function f_xi = rk4(x0, y0, x, h)

    %% Count number of iterations
    %% using step size or step
    %% height h
    n = ((x – x0) / h);

    %% Iterate for number
    %% of iterations
    y = y0;
    for i = 1:n

    %% Apply Runge Kutta
    %% Formulas to find
    %% next value of y
    k1 = h * dydx(x0, y);
    k2 = h * dydx(x0 + 0.5 *h,y + 0.5 *k1);
    k3 = h * dydx(x0 + 0.5 *h,y + 0.5 * k2);
    k4 = h * dydx(x0 + h, y + k3);

    %% Update next value of y
    f(i) = y;
    y = y + (1.0 / 6.0) * (k1 + 2*k2 + 2 *k3 + k4);
    %% Update next value of x
    x0 = x0 + h;
    f_xi(i) = y;
    end
    end
    % plot([x0:h:x],f_xi);

  • سعید خانی پرکی — says: ۳ مهر، ۱۳۹۸ در ۱:۵۰ ب٫ظ

    با سلام و تشکر بخاطر نشر مطالب علمی
    کد متلب مثال بالا رو برای استفاده دوستان میزارم :

    %% Runge Kutta method

    %% Driver method
    x0 = 0;
    y = 1;
    x = 2;
    h = 0.2;
    display(sprintf(‘The value of y at x is : %f ‘, rk4(x0, y, x, h)))

    %% A sample differential equation
    %% “dy/dx = (x – y)/2”
    function f = dydx(x,y)
    f = (x-y) / 2;
    end

    %% Finds value of y for a
    %% given x using step size h
    %% and initial value y0 at x0.
    function y = rk4(x0, y0, x, h)

    %% Count number of iterations
    %% using step size or step
    %% height h
    n = ((x – x0) / h);

    %% Iterate for number
    %% of iterations
    y = y0;
    for i = 1:n

    %% Apply Runge Kutta
    %% Formulas to find
    %% next value of y
    k1 = h * dydx(x0, y);
    k2 = h * dydx(x0 + 0.5 *h,y + 0.5 *k1);
    k3 = h * dydx(x0 + 0.5 *h,y + 0.5 * k2);
    k4 = h * dydx(x0 + h, y + k3);

    %% Update next value of y
    y = y + (1.0 / 6.0) * (k1 + 2*k2 + 2 *k3 + k4);
    %% Update next value of x
    x0 = x0 + h;
    end
    end

نظر شما چیست؟

نشانی ایمیل شما منتشر نخواهد شد.

مشاهده بیشتر