求解微分方程组的ODE算法
来源:互联网 发布:linux 平台 编辑:程序博客网 时间:2024/05/21 17:38
4阶5级龙格库塔法用于解一阶微分方程(组),对于高阶微分方程,可以将其转换为一阶微分方程组求解。原程序由John.H.Mathews编写(数值方法matlab版),但只能解微分方程,不能解微分方程组。由LiuLiu@uestc修改,使之能够解微分方程组。该程序精度比matlab自带的ode45更高。
function [Rt Rx]=rkf45(f,tspan,ya,m,tol)
% Input:
% - f function column vector
% - tspan[a,b] left & right point of [a,b]
% - ya initial value column vector
% -m initial guess for number of steps
% -tol tolerance
% Output:
% - Rt solution: vector of abscissas
% - Rx solution: vector of ordinates
% Program by John.Mathews, improved by liuliu@uestc
if length(tspan)~=2
error('length of vector tspan must be 2.');
end
if ~isnumeric(tspan)
error('TSPAN should be a vector of integration steps.');
end
if ~isnumeric(ya)
error('Ya should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.') ;
end
a=tspan(1);
b=tspan(2);
ya=ya(:);
a2 = 1/4; b2 = 1/4; a3 = 3/8; b3 = 3/32; c3 = 9/32; a4 = 12/13;
b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197; a5 = 1;
b5 = 439/216; c5 = -8; d5 = 3680/513; e5 = -845/4104; a6 = 1/2;
b6 = -8/27; c6 = 2; d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40;
r1 = 1/360; r3 = -128/4275; r4 = -2197/75240; r5 = 1/50;
r6 = 2/55; n1 = 25/216; n3 = 1408/2565; n4 = 2197/4104; n5 = -1/5;
big = 1e15;
h = (b-a)/m;
hmin = h/64;% 步长自适应范围下限
hmax = 64*h;% 步长自适应范围上限
max1 = 200;% 迭代次数上限
Y(1,:) = ya;
T(1) = a;
j = 1;
% tj = T(1);
br = b - 0.00001*abs(b);
while (T(j)<b),
if ((T(j)+h)>br), h = b - T(j); end
%caculate values of k1...k6,y1...y6
tj = T(j);
yj = Y(j,:);
y1 = yj;
k1 = h*feval(f,tj,y1);
y2 = yj+b2*k1;
if big<abs(max(y2)) return, end
k2 = h*feval(f,tj+a2*h,y2);
y3 = yj+b3*k1+c3*k2; if big<abs(max(y3)) return, end
k3 = h*feval(f,tj+a3*h,y3);
y4 = yj+b4*k1+c4*k2+d4*k3; if big<abs(max(y4)) return, end
k4 = h*feval(f,tj+a4*h,y4);
y5 = yj+b5*k1+c5.*k2+d5*k3+e5*k4; if big<abs(max(y5)) return, end
k5 = h*feval(f,tj+a5*h,y5);
y6 = yj+b6*k1+c6.*k2+d6*k3+e6*k4+f6*k5; if big<abs(max(y6)) return, end
k6 = h*feval(f,tj+a6*h,y6);
err = abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);
ynew = yj+n1*k1+n3*k3+n4*k4+n5*k5;
% error and step size control
if ( (err<tol) | (h<2*hmin) ),
Y(j+1,:) = ynew;
if ((tj+h)>br),
T(j+1) = b;
else
T(j+1) = tj + h;
end
j = j+1;
tj = T(j);
end
if (max(err)==0),
s = 0;
else
s1 = 0.84*(tol.*h./err).^(0.25);% 最佳步长值
s=min(s1);
end
if ((s<0.75)&(h>2*hmin)), h = h/2; end
if ((s>1.50)&(2*h<hmax)), h = 2*h; end
if ( (big<abs(Y(j,:))) | (max1==j) ), return, end
end
% [Rt Rx]=[T' Y];
Rt=T';
Rx=Y;
使用方法:
首先编写方程(组)文件(注意与ode45不同,这儿方程组为1Xn数组:
function dx= fun(t,x)
dx=zeros(1,2);
dx(1)=x(1)+x(2)*2+0*t;
dx(2)=3*x(1)+x(2)*2+0*t;
然后使用:
[Rt,Rx]=rkf45(@fun,[0,0.2],[6;4],100,1e-7)
- 求解微分方程组的ODE算法
- 求解一个有趣的常微分方程组
- Matlab求解微分方程组
- 全区间积分的哈明方法(常微分方程组的求解)
- 全区间积分的阿当姆斯预报校正法(常微分方程组的求解)
- 全区间积分的双边法(常微分方程组的求解)
- 全区间积分的定步长欧拉方法(常微分方程组的求解)
- Matlab求解中性类型的时滞微分方程组-中性类型的时滞微分方程
- matlab ode方程的求解
- 求解Matlab微分方程组中的时移问题!!!
- matlab求解时滞微分方程组—固定时滞
- 线性微分方程组的公式解法
- ode求解器的事件(Event)属性
- 多元一次方程组的求解
- Matlab中的PDEPE求解"瞬态型"或"发展型"非线性偏微分方程组
- matlab求解常微分方程组/传染病模型并绘制SIR曲线
- 如何用matlab绘微分方程组解的图像
- 海啸的数学模型(也是微分方程组)【作者:泰瑞陶】
- Spring基础 导入AOP+自定义标签切面
- VB 共享软件防破解设计技术初探(一)
- Linux查看CPU信息
- maven项目:java.lang.ClassNotFoundException
- 【杭电OJ】2544--最短路(最短路)
- 求解微分方程组的ODE算法
- Python 函数
- 1.5 一级指针内存模型(初学者重点)
- 红帽6.5未安装中文语言导致数据库乱码
- Ubuntu12.04(64bit)编译Android4.4源码和kernel
- 给用torch的童鞋,lua代码规范
- 学习之路 vue的介绍
- angular-cli的安装及各种坑
- 关于微信发送中文 总是显示\u***”