Lorenz方程组的MATLAB实现

首先，我们考虑Lorenz方程组

dy1/dt=sigma*(y2-y1),

dy2/dt=r*y1-y2-y1*y3,

dy3/dt=y1*y2-b*y3,

其中 sigma为普朗特数，r为瑞利数，b为非负常数，这是一个混沌系统。当我们取参数 sigma=10，r=28，sigma=8/3，初值为(0,1,0)时，采用四阶显式龙格库塔方法进行求解的MATLAB程序如下：

主程序：

clc,clear;
close all;
k = 0.01; % constant step size
tf = 100; % final time
% given initial conditions
y0(1) = 0; y0(2) = 1; y0(3) = 0;
y0 = shiftdim(y0);
% integrate
[tout,yout] = rk4(tf,k,y0);
% plot solution in phase space
figure (1)
subplot(1,3,1)
plot(yout(1,:),yout(2,:));
xlabel('y_1');
ylabel('y_2');
subplot(1,3,2)
plot(yout(1,:),yout(3,:));
xlabel('y_1');
ylabel('y_3');
subplot(1,3,3)
plot(yout(2,:),yout(3,:));
xlabel('y_2');
ylabel('y_3');
figure (2)
subplot(1,3,1)
plot(tout,yout(1,:));
xlabel('t');
ylabel('y_1');
subplot(1,3,2)
plot(tout,yout(2,:));
xlabel('t');
ylabel('y_2');
subplot(1,3,3)
plot(tout,yout(3,:));
xlabel('t');
ylabel('y_3');

子程序：显式龙格库塔法

function [tout,yout] = rk4 (tf,k,y0)
%
% [tout,yout] = rk4 (tf,k,y0)
% solves non-stiff IVODEs based on classical 4-stage Runge-Kutta
%
% y’ = ffun (t,y), y(0) = y0, 0 < t < tf
% ffun returns a vector of m components for a given
% t and y of m components.
%
% Input:
% tf - final time
% k - uniform step size
% y0 - initial values y(0) (size (m,1))
%
% Output:
% tout - times where solution is computed
% yout - solution vector at times tout.
% prepare
N = tf / k;
tout = [0];
fo = ffun (tout(1),y0);
yout(:,1) = y0;
% loop in time
for n=1:N
t = tout(n);
y = yout(:,n);
Y1 = y;
K1 = ffun(t,Y1);
Y2 = y + k/2*K1;
t = t + k/2;
K2 = ffun(t,Y2);
Y3 = y + k/2*K2;
K3 = ffun(t,Y3);
Y4 = y + k*K3;
t = t + k/2;
K4 = ffun(t,Y4);
tout(n+1) = t;
yout(:,n+1) = y + k/6*(K1+2*K2+2*K3+K4);
end
end

子程序：求解函数

function f = ffun(t,y)
% function defining ODE
sigma = 10; b = 8/3; r = 28;
f = zeros(3,1);
f(1) = sigma *(y(2)-y(1));
f(2) = r * y(1) - y(2) - y(1)*y(3);
f(3) = y(1)*y(2) - b*y(3);

结果：

 

参考文献：

[1] Ascher, Uri M . Numerical Methods for Evolutionary Differential Equations, SIAM, 2008.