标签:com 名称 style end isp floor 算法 height 维数
非刚性常微分方程的数值解法通常会用四阶龙格库塔算法,其matlab函数对应ode45。
对于dy/dx = f(x,y),y(0)=y0。
其四阶龙格库塔公式如下:
对于通常计算,四阶已经够用,四阶以上函数f(x,y)计算工作量大大增加而精度提高较慢。
下面以龙格库塔法解洛伦兹方程为例:
matlab代码如下:
main.m:
clear all; close all; clc; %系统龙格库塔法 [t,h] = ode45(@test_fun,[0 40],[12 4 0]); plot3(h(:,1),h(:,2),h(:,3)); grid on; %自定义龙格库塔法 [t1,h1]=runge_kutta(@test_fun,[12 4 0],0.01,0,40); figure; plot3(h1(1,:),h1(2,:),h1(3,:),‘r‘) grid on;
runge_kutta.m(函数参考网络):
%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数) function [x,y]=runge_kutta(ufunc,y0,h,a,b) n=floor((b-a)/h); %步数 x(1)=a; %时间起点 y(:,1)=y0; %赋初值,可以是向量,但是要注意维数 for i=1:n %龙格库塔方法进行数值求解 x(i+1)=x(i)+h; k1=ufunc(x(i),y(:,i)); k2=ufunc(x(i)+h/2,y(:,i)+h*k1/2); k3=ufunc(x(i)+h/2,y(:,i)+h*k2/2); k4=ufunc(x(i)+h,y(:,i)+h*k3); y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6; end
test_fun(洛伦兹方程):
%构造微分方程 function dy=test_fun(t,y) a = 16; b = 4; c = 45; dy=[a*(y(2)-y(1)); c*y(1)-y(1)*y(3)-y(2); y(1)*y(2)-b*y(3)];
得到很经典的洛伦兹吸引子,结果如下:
爸妈这电脑只有Matlab2010,画出来的图好丑。。。
参考:
https://wenku.baidu.com/view/8211fbd428ea81c758f57893.html
标签:com 名称 style end isp floor 算法 height 维数
原文地址:https://www.cnblogs.com/tiandsp/p/12238024.html