matlab：使用四阶龙格库塔方法求解微分方程组

%书籍：常用数值算法及其matlab实现
%第10章 常微分方程初值问题的数值解法,例10.14使用
%四阶龙格库塔方法
function [t,z] = rk4symeq(fun, t0, tf, Za, h)
%fun:微分方程的右表达式
%t0, tn为区间
%Za为初值,是列向量
M = floor(tf-t0)/h ;      %离散点的个数M+1
if t0 >= tf
    printf('左端点必须小于右端点');
    return;
end
N = length(Za);           %获得变量个数,N
z = zeros(M+1, N);
t = zeros(M +1, 1);
t =[t0 : h :tf]';
z(1,:) = Za';            %假设Za为列向量，与微分方程中的变量方向统一，变成行向量

for i = 1:M
    K1 =  feval(fun, t(i) , z(i,:));                    %K是行向量
    K2 =  feval(fun, t(i)+1/2*h ,z(i,:)+1/2* h*K1);
    K3 =  feval(fun, t(i)+1/2*h ,z(i,:)+1/2* h*K2);
    K4 =  feval(fun, t(i)+ h ,z(i,:)+ h*K3);   
    z(i+1,:) = z(i,:) +h/6 *(K1 + 2*K2 + 2*K3 + K4);
end
方程组1-输出为行向量
%书籍：常用数值算法及其matlab实现
%第10章 常微分方程初值问题的数值解法
%四阶龙格库塔，例10.15
function s = exa10_15(t,z)
%z是个向量,1*3
%输出s也是向量,1*3
s = zeros(1,3);
dy1 = -10*z(1) +10*z(2);
dy2 = 28*z(1) - z(2) - z(1)*z(3);
dy3 = -8/3*z(3) + z(1)*z(2);
s = [dy1 dy2 dy3];          %输出s为行向量

微分方程组——输出为列向量
%书籍：常用数值算法及其matlab实现
%第10章 常微分方程初值问题的数值解法
%四阶龙格库塔，例10.15
function s = exa10_15b(t,z)
%z是个向量,1*3
%输出s也是向量,1*3
s = zeros(1,3);
dy1 = -10*z(1) +10*z(2);
dy2 = 28*z(1) - z(2) - z(1)*z(3);
dy3 = -8/3*z(3) + z(1)*z(2);
s = [dy1; dy2; dy3];          %输出s为列向量，用于ode45调用
主函数
clear all;clc;close all;
%书籍：常用数值算法及其matlab实现
%第10章 常微分方程组初值问题的数值解法
%四阶龙格库塔，例10.15
%函数原型 function [t,z] = rk4symeq(fun, t0, tn, Za, h)
format long
t0 = 0; tf = 5;
Za = [-8 ; 8; 27];     %x初值
h1 = 0.01; 

[t1,z1] = rk4symeq(@exa10_15, t0, tf , Za, h1);
figure(1)
plot(t1,z1(:,1),'b',t1,z1(:,2), 'r',t1,z1(:,3), 'g--')
legend('y1','y2','y3')
figure(2)
plot3(z1(:,1),z1(:,2),z1(:,3));
xlabel('x');ylabel('y');zlabel('z');

[t2,z2] = ode45(@exa10_15b, [0,5] ,[-8 , 8 ,27]);        %用系统函数
figure(3)
plot(t2,z2(:,1),'b',t2,z2(:,2), 'r',t2,z2(:,3), 'g--')
legend('y1','y2','y3')

运行结果：