نحوه حل ODE ها در MATLAB با استفاده از Runge-Kutta

آموزش کدنویسی الگوریتم مرتبه چهارم Runge-Kutta برای حل معادلات دیفرانسیل معمولی

معادلات دیفرانسیل معمولی (ODE) در تمام جنبه های مهندسی و علوم رایج است. آنها برای مدل‌سازی پدیده‌های فیزیکی مانند جریان هوای متلاطم، مدارهای ماهواره‌ای، انتقال حرارت، ارتعاشات و موارد بی‌شماری استفاده می‌شوند.

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

معادلات دیفرانسیل معمولی

در این مقاله، فرض می‌کنم که شما درک اساسی از حساب دیفرانسیل و انتگرال و معادلات دیفرانسیل دارید. اصطلاح معادلات دیفرانسیل معمولی به معادلات دیفرانسیل اطلاق می شود که شامل توابع یک متغیر مستقل و مشتقات آن توابع است. برای اکثر ODE های مهندسی و علوم، متغیر مستقل زمان است و معادلات به شکل زیر است:

انجام پروژه متلب در فریلنس پروژه

در اینجا، بردار y می تواند حالت فضایی، جرم، دما و غیره باشد و t زمان است. متأسفانه، اکثر ODE ها را نمی توان به صورت تحلیلی حل کرد. این بدان معنی است که نمی توان تابعی را برای ارائه یک راه حل دقیق برای معادله دیفرانسیل استخراج کرد. راه دیگر برای حل این معادلات دیفرانسیل استفاده از روش های انتگرال گیری عددی است.

ادغام عددی

اصطلاح ادغام عددی برای اولین بار در سال 1915 ابداع شد، اما مزایای آن تا زمانی که کامپیوترهای مدرن به طور واقعی مشاهده نشدند. انتگرال گیری عددی روشی برای تقریب تغییر تابع y در طول زمان با دانستن معادلات دیفرانسیل حاکم بر تغییر y در زمان است. آنها همانطور که گفته شد یک تخمین هستند، بنابراین فقط به اندازه روش ها و تکنیک های مورد استفاده خوب هستند. یکی از رایج‌ترین شکل‌های ادغام عددی، روش Runge-Kutta مرتبه چهارم (یا RK4) است.

RK4 را می توان با معادلات و نمودار زیر توصیف کرد. بیشتر اوقات، شما یک فرم برداری از ODE ها برای حل خواهید داشت، بنابراین شکل برداری RK4 نشان داده می شود. در معادلات، مقادیر k تخمین شیب y هستند که با استفاده از معادلات دیفرانسیل در مکان های نشان داده شده در نمودار محاسبه می شوند. هنگامی که مقادیر k تعیین می شوند، میانگین می شوند و در گام زمانی، h ضرب می شوند. سپس این مقدار به مقدار فعلی y اضافه می شود تا تقریب مرحله زمانی بعدی y بدست آید. این فرآیند تا رسیدن به تعداد مراحل زمانی مورد نظر تکرار می شود. توجه داشته باشید که هر چه مقدار h استفاده شده کوچکتر باشد، نتیجه به مقدار واقعی دقیق تر خواهد بود.

اسکریپت MATLAB RK4

برای این مثال، فرض کنید یک ذره داریم که حرکت آن با معادلات دیفرانسیل معمولی زیر توصیف می شود:

این مجموعه از ODE ها راحت است زیرا می توان آن را به صورت تحلیلی حل کرد (راه حل نشان داده شده در زیر). اینها انتخاب شدند تا الگوریتم RK4 ما بتواند با راه حل واقعی مقایسه شود. به خاطر داشته باشید که هنگام استفاده از ODE های خود، آنها ممکن است به صورت تحلیلی قابل حل نباشند، بنابراین ممکن است چیزی برای مقایسه نتایج خود نداشته باشید.

می‌توانیم کد را با تعریف یک تابع، مدل شروع کنیم، که وقتی مقادیر t، v و w به آن ارسال می‌شود، مقدار ODE‌ها را برمی‌گرداند. مطمئن شوید که این تعریف تابع را در پایین اسکریپت خود قرار داده اید.

function dydt = model(t,y)
    v = y(1);
    w = y(2);
    dv_dt = w;
    dw_dt = 6*v-w;
    dydt = [dv_dt,dw_dt];
end

بیایید اسکریپت اصلی را اکنون شروع کنیم. من همیشه آن را با پاک کردن پنجره فرمان، پاک کردن متغیرها و بستن ارقام شروع می کنم.

clc
clear variables
close all

برای استفاده از روش RK4، باید شرایط اولیه، y0، و یک آرایه زمانی، t را تعریف کنیم. در اینجا، t تعریف شده است که از 0 به 1 در 1000 افزایش می رود و شرط اولیه برای v 3 و w 1 است. آرایه زمانی برای ایجاد h و بعد از آن، مقادیر دقیق برای v و w استفاده می شود.

% Initial Conditions and Simulation Time
y0 = [3, 1];  % y0 = [v0, w0]
t = linspace(0,1,1000)';

سپس، می‌توانیم یک آرایه، yₖᵤₜₜₐ ایجاد کنیم تا تقریب‌های RK4 را در هر مرحله زمانی ذخیره کنیم. با استفاده از شرایط اولیه برای v و w، می توانیم اولین اندیس (i = 1، t = 0) آرایه RK4 را تعریف کنیم. همچنین می‌توانیم متغیر مرحله زمانی، h را تعریف کنیم. سپس یک حلقه for برای تکرار در طول آرایه زمانی معرفی می کنیم. در هر تکرار یک مقدار جدید به yₖᵤₜₜₐ اضافه می کنیم و سپس از این تکرار برای ایجاد مقدار دیگری برای yₖᵤₜₜₐ استفاده می کنیم. این روند تا زمانی که به آخرین مرحله زمانی برسیم (t = 1) تکرار خواهد شد.

% Runge-Kutta 4th-Order Algorithm
y_Kutta = zeros(length(t), 2);
y_Kutta(1, 🙂 = y0;
h = t(2)-t(1);  % Constant time step
for i = 2:length(t)
    k1 = model(t(i-1), y_Kutta(i-1, :));
    k2 = model(t(i-1)+h/2, y_Kutta(i-1, :)+k1*h/2);
    k3 = model(t(i-1)+h/2, y_Kutta(i-1, :)+k2*h/2);
    k4 = model(t(i-1)+h, y_Kutta(i-1, :)+k3*h);
    y_Kutta(i, 🙂 = y_Kutta(i-1, :)+(k1/6+k2/3+k3/3+k4/6)*h;
end

همانطور که قبلا ذکر شد، این مجموعه از ODE ها یک راه حل تحلیلی دارد، بنابراین بیایید یک آرایه از آن مقادیر ایجاد کنیم و آن را با نتیجه خود از RK4 مقایسه کنیم. برای مقایسه

e، ما به سادگی می توانیم تفاوت بین راه حل دقیق و جواب تقریبی خود را محاسبه کنیم.

% Exact
y_Exact = 2*exp(2*t).*[1, 2]-exp(-3.*t)*[-1, 3];% Calculating the Difference Between Exact and RK4 Solutions
diff_Kutta = y_Exact-y_Kutta;

در نهایت، ما تمام داده هایی را که برای مقایسه روش مرتبه چهارم Runge-Kutta با راه حل دقیق نیاز داریم، در اختیار داریم. مرحله آخر ترسیم نتایج است.

% Comparison Subplot
figure; subplot(2,1,1);  hold on;
plot(y_Exact(:, 1), y_Exact(:, 2))
plot(y_Kutta(:, 1), y_Kutta(:, 2))
title('RK4 Compared to the Exact Solution', 'Interpreter', 'Latex')
xlabel('v', 'Interpreter', 'Latex')
ylabel('w', 'Interpreter', 'Latex')
legend('Exact', 'RK4', 'Location', 'Best')
grid on; hold off% Difference Subplot
subplot(2,1,2); hold on;
plot(t, diff_Kutta(:, 1))
plot(t, diff_Kutta(:, 2))
xlabel('t', 'Interpreter', 'Latex')
ylabel('Difference', 'Interpreter', 'Latex')
legend('RK4 [v]', 'RK4 [w]', 'Location', 'Best')
grid on; hold off

نتایج:

نتایج Range-Kutta مرتبه چهارم [ایجاد شده توسط نویسنده]

همانطور که می بینید، نتایج روش Runge-Kutta مرتبه 4 برای این مثال بسیار دقیق است. نتایج RK4 از جواب دقیق در نمودار بالا قابل تشخیص نیستند و تفاوت بین نتایج برای v و w و جواب دقیق عملاً ناچیز است. Runge-Kutta را می توان فراتر از مرتبه 4 برای کاهش خطا گسترش داد، اما این احتمالاً برای اکثر برنامه ها ضروری نیست.

با تشکر از شما برای خواندن مقاله! امیدوارم کمی در مورد RK4 یاد گرفته باشید! من را دنبال کنید و مقالات دیگر من در مورد پایتون، ریاضیات و مکانیک مداری را بررسی کنید! اگر نظر یا پیشنهادی دارید، به من بگویید!
نحوه حل ODE ها در MATLAB