4 详解matlab实现龙格库塔算法求解复杂常微分方程组

目录

4.1 题目

4.2 问题背景

4.3 算法

4.4 matlab编程实现

4.5 结果展示

4.6 讨论

---

4.1 题目

4.2 问题背景

 

4.3 算法

 

4.4 matlab编程实现

主程序：

clear;

fun1=@(x,y1,y2,y3)-0.013*y1-1000*y1*y2;  % y1(x)的一阶导数d(y1(x))/dx

fun2=@(x,y1,y2,y3)-2500*y2*y3; 

fun3=@(x,y1,y2,y3)-0.013*y1-1000*y1*y2-2500*y2*y3;

x0=0;y0=[1;1;0];h=0.0001;       %赋初值

xn1=0.001;xn2=0.004;xn3=0.02;   %在xn1处的待求函数值y(xn1)

yn1=RungeK4(x0,xn1,y0,h,fun1,fun2,fun3)

yn2=RungeK4(x0,xn2,y0,h,fun1,fun2,fun3)

yn3=RungeK4(x0,xn3,y0,h,fun1,fun2,fun3)

函数：

function [yn,Y]=RungeK(x0,xn,y0,h,fun1,fun2,fun3)

%四阶古典龙格库塔算法求解1元-3元的微分方程通用解法

% x0为初始起点， x1为待求点自变量， h为步长， y为待求点函数值 , yn为所求点的三个函数值组成的向量,Y为中间值矩阵

format long

m=length(y0);  %m表示y向量包含的函数个数

n=(xn-x0)/h+1;

Y=zeros(3,n);   %创建二维数组以容纳n个列向量y(x,i)=[y1(x,i);y2(x,i);y3(x,i)]

yn=zeros(m,1);

Y(1:m,1)=y0;  %赋初值

K1=zeros(3,n);K2=zeros(3,n);K3=zeros(3,n);K4=zeros(3,n);   %创建二维数组以容纳n个列向量,K1（1，i）对应y1的k1(i)

for i=1:n

    x(i)=x0+(i-1)*h;

   

    %为向量K1,K2,K3,K4依次赋值

    K1(1,i)=fun1(x(i),Y(1,i),Y(2,i),Y(3,i));  %其对应y1的k1（i）

    K1(2,i)=fun2(x(i),Y(1,i),Y(2,i),Y(3,i));  %其对应y2的k1（i）

    K1(3,i)=fun3(x(i),Y(1,i),Y(2,i),Y(3,i));  %其对应y3的k1（i）

   

    K2(1,i)=fun1(x(i)+1/2*h,Y(1,i)+1/2*K1(1,i)*h,Y(2,i)+1/2*K1(1,i)*h,Y(3,i)+1/2*K1(1,i)*h);  %为向量K2赋值，其对应y1的k2（i）

K2(2,i)=fun2(x(i)+1/2*h,Y(1,i)+1/2*K1(2,i)*h,Y(2,i)+1/2*K1(2,i)*h,Y(3,i)+1/2*K1(2,i)*h);

    K2(3,i)=fun3(x(i)+1/2*h,Y(1,i)+1/2*K1(3,i)*h,Y(2,i)+1/2*K1(3,i)*h,Y(3,i)+1/2*K1(3,i)*h);

   

    K3(1,i)=fun1(x(i)+1/2*h,Y(1,i)+1/2*K2(1,i)*h,Y(2,i)+1/2*K2(1,i)*h,Y(3,i)+1/2*K2(1,i)*h);  %为向量K3赋值，其对应y1的k3（i）

    K3(2,i)=fun2(x(i)+1/2*h,Y(1,i)+1/2*K2(2,i)*h,Y(2,i)+1/2*K2(2,i)*h,Y(3,i)+1/2*K2(2,i)*h);

    K3(3,i)=fun3(x(i)+1/2*h,Y(1,i)+1/2*K2(3,i)*h,Y(2,i)+1/2*K2(3,i)*h,Y(3,i)+1/2*K2(3,i)*h);

   

    K4(1,i)=fun1(x(i)+h,Y(1,i)+K3(1,i)*h,Y(2,i)+K3(1,i)*h,Y(3,i)+K3(1,i)*h);  %为向量K4赋值，其对应y1的k4（i）

    K4(2,i)=fun2(x(i)+h,Y(1,i)+K3(2,i)*h,Y(2,i)+K3(2,i)*h,Y(3,i)+K3(2,i)*h);

    K4(3,i)=fun3(x(i)+h,Y(1,i)+K3(3,i)*h,Y(2,i)+K3(3,i)*h,Y(3,i)+K3(3,i)*h);

   

    for j=1:3

        %为向量y(x)依次赋值

        Y(j,i+1)=Y(j,i)+h/6*(K1(j,i)+2*K2(j,i)+2*K3(j,i)+K4(j,i));

    end

%     Y(1,i+1)=Y(1,i)+h/6*(K1(1,i)+2*K2(1,i)+2*K3(1,i)+K4(1,i));

%     Y(2,i+1)=Y(2,i)+h/6*(K1(2,i)+2*K2(2,i)+2*K3(2,i)+K4(2,i));

%     Y(3,i+1)=Y(3,i)+h/6*(K1(3,i)+2*K2(3,i)+2*K3(3,i)+K4(3,i));

end

for j=1:m

    %为所求函数值向量yn依次赋值

    yn(j)=Y(j,n);

end

   

   

4.5 结果展示

>> work4

yn1 =

   0.329092932578304

   1.442413401551611

  -0.189667611588121

yn2 =

   0.001386837752709

   1.943934950799432

  -0.001112694535328

yn3 =

   0.000000000000000

   1.946359874539511

  -0.000000000000000

4.6 讨论

（1）通过对y(x)、k1、k2、k3、k4向量化可以同时计算多元一阶微分方程，较为简便

（2）龙格库塔法比欧拉法精度较高，且收敛快，充分光滑时阶数越高精度越高。

（3）由于龙格库塔法基于Taylor公式，因此它对光滑性有一定要求，若光滑性差，则四阶龙格库塔法数值解精度可能不如二阶龙格库塔法。

（4）通过课本例题验证，本程序通用性较强且误差很小