龍格庫塔法是一種求解高階常微分方程的常用方法，在工程當(dāng)中應(yīng)用廣泛，例如求解物體的運動方程等。
這里我們通過matlab程序編寫龍格庫塔算法求解二元常微分方程組，假設(shè)有常微分方程組：
$?\begin{array}{l}\ddot{x}?\dot{x}+2\ddot{y}+\dot{y}=?2\mathrm{s}\mathrm{i}\mathrm{n}t?3\mathrm{c}\mathrm{o}\mathrm{s}t\\ \ddot{x}+\ddot{y}=?\mathrm{s}\mathrm{i}\mathrm{n}t?\mathrm{c}\mathrm{o}\mathrm{s}t\\ x(0)=0\\ y(0)=1\\ \dot{x}(0)=1\\ \dot{y}(0)=0\end{array}$
該方程有精確解：
$\{\begin{array}{l}x=\mathrm{s}\mathrm{i}\mathrm{n}t\\ y=\mathrm{c}\mathrm{o}\mathrm{s}t\end{array}$
令$u_1=x$,$u_2=\dot{x}$,$w_1=y$,$w_2=\dot{y}$,則原式可寫為
$?\begin{array}{l}f1(t)=u_2=\dot{x}\\ f2(t)=u_2^{\prime }=\ddot{x}=\mathrm{c}\mathrm{o}\mathrm{s}t?\dot{x}+\dot{y}\\ f3(t)=w_2=\dot{y}\\ f4(t)=w_2^{\prime }=\ddot{y}=?\mathrm{s}\mathrm{i}\mathrm{n}t?2\mathrm{c}\mathrm{o}\mathrm{s}t+\dot{x}?\dot{y}\end{array}$
可以編寫四個子函數(shù)如下：

function output = f1(x,u1,u2,w1,w2) output = u2; function output = f2(x,u1,u2,w1,w2) output = cos(x)-u2+w2; function output = f3(x,u1,u2,w1,w2) output = w2; function output = f4(x,u1,u2,w1,w2) output = -sin(x)-2*cos(x)-w2+u2;

四階龍格庫塔函數(shù)程序

function [u1,u2,w1,w2] = RK4_2variable(u1,u2,w1,w2,h,a,b)x = a:h:b;for i = 1:length(x)-1k11 = f1(x(i) , u1(i) , u2(i) , w1(i) , w2(i)); k21 = f2(x(i) , u1(i) , u2(i) , w1(i) , w2(i)); L11 = f3(x(i) , u1(i) , u2(i) , w1(i) , w2(i)); L21 = f4(x(i) , u1(i) , u2(i) , w1(i) , w2(i));k12 = f1(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2); k22 = f2(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2); L12 = f3(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2); L22 = f4(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);k13 = f1(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2); k23 = f2(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2); L13 = f3(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2); L23 = f4(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);k14 = f1(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23); k24 = f2(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23); L14 = f3(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23); L24 = f4(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);u1(i+1) = u1(i) + h/6 * (k11 + 2*k12 + 2*k13 + k14); u2(i+1) = u2(i) + h/6 * (k21 + 2*k22 + 2*k23 + k24); w1(i+1) = w1(i) + h/6 * (L11 + 2*L12 + 2*L13 + L14); w2(i+1) = w2(i) + h/6 * (L21 + 2*L22 + 2*L23 + L24);end end

計算主程序

% main func clear;clc; u1(1) = 0; u2(1) = 1; w1(1) = 1; w2(1) = 0; h=0.01; a = 0;b=20; [u1,u2,w1,w2] = RK4_2variable(u1,u2,w1,w2,h,a,b); figure plot(a:h:b,u1,'r-'); hold on plot(a:h:b,sin(a:h:b),'b-.'); xlabel('time'); ylabel('x'); legend('計算值','精確值');

計算結(jié)果如圖

總結(jié)

以上是生活随笔為你收集整理的Matlab 四阶龙格库塔法求解二元常微分方程组的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。