㈠ 龍格庫塔法求解微分方程,matlab怎麼編程
龍格庫塔法是一種高效的數值方法,用於求解常微分方程。這里提供4階龍格庫塔法的MATLAB代碼示例:
function [Y] = RK45(t,X,f,h)
K1 = f(t,X);
K2 = f(t+h/2,X+h/2*K1);
K3 = f(t+h/2,X+h/2*K2);
K4 = f(t+h,X+h*K3);
Y = X + h/6*(K1 + 2*K2 + 2*K3 + K4);
該函數接受當前時間t、狀態向量X、微分方程f和步長h作為輸入參數。通過計算四個斜率值K1、K2、K3、K4,進而得到新的狀態向量Y。
接下來定義一個具體的微分方程函數f,例如:
function dxdt = f(t,x)
dxdt(1) = exp(x(1)*sin(t)) + x(2);
dxdt(2) = exp(x(2)*cos(t)) + x(1);
dxdt = dxdt(:);
這個函數描述了兩個狀態變數x1和x2隨時間t的變化關系。x(1)對應x1,x(2)對應x2。這里採用逐項賦值的方式,最後將結果轉化為列向量。
在確定了初始條件和參數後,可以使用以下代碼求解t0到t1的時間區間內的軌跡:
t0 = 0; t1 = 5; h = 0.02; x0 = [-1; -1];
T = t0:h:t1;
X = zeros(length(x0), length(T));
X(:, 1) = x0;
for j = 1:length(T)-1
X(:, j+1) = RK45(T(j), X(:, j), @(t, x) f(t, x), h);
end
plot(T, X(1, :));
hold on;
plot(T, X(2, :), 'r');
這段代碼中,首先定義了初始時間t0、終止時間t1、步長h和初始狀態向量x0。然後通過循環調用RK45函數,逐步更新狀態向量X,最終得到時間序列T和狀態序列X。最後,利用plot函數繪制了兩個狀態變數隨時間的變化曲線。