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程序如下:
主程序:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | 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' ); |
子程序:显式龙格库塔法
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 | 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 |
子程序:求解函数
1 2 3 4 5 6 7 | 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.
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· Java 中堆内存和栈内存上的数据分布和特点
· 开发中对象命名的一点思考
· .NET Core内存结构体系(Windows环境)底层原理浅谈
· C# 深度学习:对抗生成网络(GAN)训练头像生成模型
· .NET 适配 HarmonyOS 进展
· 用 DeepSeek 给对象做个网站,她一定感动坏了
· DeepSeek+PageAssist实现本地大模型联网
· 手把手教你更优雅的享受 DeepSeek
· 腾讯元宝接入 DeepSeek R1 模型,支持深度思考 + 联网搜索,好用不卡机!
· 从 14 秒到 1 秒:MySQL DDL 性能优化实战