✅作者简介热爱科研的Matlab仿真开发者擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。完整代码获取 定制创新 论文复现点击Matlab科研工作室 关注我领取海量matlab电子书和数学建模资料个人信条做科研博学之、审问之、慎思之、明辨之、笃行之是为博学慎思明辨笃行。 内容介绍一、引言在科学与工程领域混沌系统广泛存在从气象学中的天气预测到电子电路中的信号处理都能发现其身影。传统整数阶混沌系统已得到深入研究但分数阶混沌系统由于其记忆性和遗传性能更精确地描述许多复杂现象。基于数据驱动的分数阶混沌系统建模通过对实际观测数据的分析构建分数阶混沌模型为理解和预测复杂动态系统提供了新的视角和方法。二、分数阶混沌系统基础二分数阶混沌系统特点与整数阶混沌系统相比分数阶混沌系统具有独特的动力学特性。由于分数阶导数的记忆性分数阶混沌系统对初始条件的敏感性更强其相空间轨迹更加复杂。例如在一些分数阶混沌电路中通过调整分数阶的阶数可以观察到不同程度的混沌行为且系统的吸引子形态也会发生变化。这种灵活性为模拟复杂的自然现象和工程系统提供了更强大的工具。三、数据驱动建模的必要性与优势一传统建模方法的局限传统的混沌系统建模通常基于物理原理和数学推导需要对系统的内部结构和参数有深入了解。然而在许多实际场景中系统的物理机制可能非常复杂难以精确描述或者系统参数难以直接测量。例如在生物系统中由于生物过程的高度复杂性和不确定性基于物理原理的建模方法往往无法准确刻画系统的动态行为。二数据驱动建模优势无需精确物理模型数据驱动建模方法直接从观测数据出发无需预先知道系统的物理模型。通过对大量数据的分析和挖掘能够自动发现数据中的潜在规律构建系统模型。这使得数据驱动建模适用于各种复杂系统尤其是那些物理机制不明确的系统。适应复杂动态变化实际系统往往处于动态变化中数据驱动建模可以根据新的数据不断更新模型适应系统的变化。例如在电力系统中负荷需求会随着时间、季节等因素变化基于数据驱动的建模方法能够实时调整模型准确预测电力负荷的变化。四、基于数据驱动的分数阶混沌系统建模步骤一数据采集确定数据源根据研究对象确定合适的数据采集方式和数据源。对于物理实验系统可以通过传感器实时采集系统的输出数据如在混沌电路实验中采集电路的电压、电流等信号。对于一些实际应用场景如气象数据可以从气象站、卫星观测等渠道获取历史数据。保证数据质量采集到的数据应具有足够的长度和精度以反映系统的动态特性。同时要对数据进行预处理去除噪声和异常值。可以采用滤波算法如卡尔曼滤波对数据进行平滑处理提高数据的质量。二特征提取混沌特征分析利用混沌理论中的方法如相空间重构、李雅普诺夫指数计算等分析采集数据的混沌特性。相空间重构通过延迟坐标法将一维时间序列数据映射到高维空间恢复系统的动力学信息。李雅普诺夫指数则用于衡量系统对初始条件的敏感性正的李雅普诺夫指数是混沌系统的重要特征之一。分数阶特征挖掘从数据中挖掘与分数阶相关的特征。由于分数阶系统具有记忆性数据的自相关性在分数阶特征提取中具有重要作用。可以通过分析数据的自相关函数在不同时间尺度上的衰减特性提取与分数阶相关的信息。三模型构建选择建模方法常见的数据驱动建模方法包括神经网络、支持向量机回归等。神经网络具有强大的非线性拟合能力能够逼近任意复杂的函数关系。支持向量机回归则在小样本数据建模中表现出色能够有效避免过拟合问题。根据数据特点和建模需求选择合适的建模方法。确定模型结构以神经网络为例需要确定网络的层数、每层的神经元个数等结构参数。可以通过交叉验证、网格搜索等方法优化模型结构提高模型的泛化能力。对于分数阶混沌系统建模还需要考虑如何将分数阶微积分引入模型中例如通过在神经网络的隐藏层中加入分数阶微分算子构建分数阶神经网络模型。四模型验证与优化模型验证将采集的数据划分为训练集、验证集和测试集。使用训练集对模型进行训练利用验证集调整模型的超参数最后在测试集上验证模型的性能。常用的验证指标包括均方误差MSE、平均绝对误差MAE等这些指标用于衡量模型预测值与实际值之间的偏差。模型优化根据验证结果对模型进行优化。如果模型在测试集上的误差较大可以尝试调整模型结构、增加数据量或采用更复杂的特征提取方法。同时也可以结合其他优化算法如遗传算法、粒子群优化算法等对模型的参数进行全局优化提高模型的性能。⛳️ 运行结果 部分代码function [t, y] fde12(alpha,fdefun,t0,tfinal,y0,h,param,mu,mu_tol)%FDE12 Solves an initial value problem for a non-linear differential% equation of fractional order (FDE). The code implements the% predictor-corrector PECE method of Adams-Bashforth-Moulton type% described in [1].%% [T,Y] FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,h) integrates the initial value% problem for the FDE, or the system of FDEs, of order ALPHA 0% D^ALPHA Y(t) FDEFUN(T,Y(T))% Y^(k)(T0) Y0(:,k1), k0,...,m-1% where m is the smallest integer grater than ALPHA and D^ALPHA is the% fractional derivative according to the Caputos definition. FDEFUN is a% function handle corresponding to the vector field of the FDE and for a% scalar T and a vector Y, FDEFUN(T,Y) must return a column vector. The% set of initial conditions Y0 is a matrix with a number of rows equal to% the size of the problem (hence equal to the number of rows of the% output of FDEFUN) and a number of columns depending on ALPHA and given% by m. The step-size H0 is assumed constant throughout the integration.%% [T,Y] FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM) solves as above with% the additional set of parameters for the FDEFUN as FDEFUN(T,Y,PARAM).%% [T,Y] FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU) solves the FDE with% the selected number MU of multiple corrector iterations. The following% values for MU are admissible:% MU 0 : the corrector is not evaluated and the solution is provided% just by the predictor method (the first order rectangular rule);% MU 0 : the corrector is evaluated by the selected number MU of times;% the classical PECE method is obtained for MU1;% MU Inf : the corrector is evaluated for a certain number of times% until convergence of the iterations is reached (for convergence the% difference between two consecutive iterates is tested).% The defalut value for MU is 1%% [T,Y] FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU,MU_TOL) allows to% specify the tolerance for testing convergence when MU Inf. If not% specified, the default value MU_TOL 1.E-6 is used.%% FDE12 is an implementation of the predictor-corrector method of% Adams-Bashforth-Moulton studied in [1]. Convergence and accuracy of% the method are studied in [2]. The implementation with multiple% corrector iterations has been proposed and discussed for multiterm FDEs% in [3]. In this implementation the discrete convolutions are evaluated% by means of the FFT algorithm described in [4] allowing to keep the% computational cost proportional to N*log(N)^2 instead of N^2 as in the% classical implementation; N is the number of time-point in which the% solution is evaluated, i.e. N (TFINAL-T)/H. The stability properties% of the method implemented by FDE12 have been studied in [5].%% [1] K. Diethelm, A.D. Freed, The Frac PECE subroutine for the numerical% solution of differential equations of fractional order, in: S. Heinzel,% T. Plesser (Eds.), Forschung und Wissenschaftliches Rechnen 1998,% Gessellschaft fur Wissenschaftliche Datenverarbeitung, Gottingen, 1999,% pp. 57-71.%% [2] K. Diethelm, N.J. Ford, A.D. Freed, Detailed error analysis for a% fractional Adams method, Numer. Algorithms 36 (1) (2004) 31-52.%% [3] K. Diethelm, Efficient solution of multi-term fractional% differential equations using P(EC)mE methods, Computing 71 (2003), pp.% 305-319.%% [4] E. Hairer, C. Lubich, M. Schlichte, Fast numerical solution of% nonlinear Volterra convolution equations, SIAM J. Sci. Statist. Comput.% 6 (3) (1985) 532-541.%% [5] R. Garrappa, On linear stability of predictor-corrector algorithms% for fractional differential equations, Internat. J. Comput. Math. 87% (10) (2010) 2281-2290.%% Copyright (c) 2011-2012, Roberto Garrappa, University of Bari, Italy% garrappa at dm dot uniba dot it% Revision: 1.2 - Date: July, 6 2012% Check inputsif nargin 9mu_tol 1.0e-6 ;if nargin 8mu 1 ;if nargin 7param [] ;endendend% Check order of the FDEif alpha 0error(MATLAB:fde12:NegativeOrder,...[The order ALPHA of the FDE must be positive. The value ...ALPHA %f can not be accepted. See FDE12.], alpha);end% Check the step--size of the methodif h 0error(MATLAB:fde12:NegativeStepSize,...[The step-size H for the method must be positive. The value ...H %e can not be accepted. See FDE12.], h);end% Structure for storing initial conditionsic.t0 t0 ;ic.y0 y0 ;ic.m_alpha ceil(alpha) ;ic.m_alpha_factorial factorial(0:ic.m_alpha-1) ;% Structure for storing information on the problemProbl.ic ic ;Probl.fdefun fdefun ;Probl.problem_size size(y0,1) ;Probl.param param ;% Check number of initial conditionsif size(y0,2) ic.m_alphaerror(MATLAB:fde12:NotEnoughInputs, ...[A not sufficient number of assigned initial conditions. ...Order ALPHA %f requires %d initial conditions. See FDE12.], ...alpha,ic.m_alpha);end% Check compatibility size of the problem with size of the vector fieldf_temp f_vectorfield(t0,y0(:,1),Probl) ;if Probl.problem_size ~ size(f_temp,1)error(MATLAB:fde12:SizeNotCompatible, ...[Size %d of the problem as obtained from initial conditions ...(i.e. the number of rows of Y0) not compatible with the ...size %d of the output of the vector field FDEFUN. ...See FDE12.], Probl.problem_size,size(f_temp,1));end% Number of points in which to evaluate weights and solutionr 16 ;N ceil((tfinal-t0)/h) ;Nr ceil((N1)/r)*r ;Q ceil(log2(Nr/r)) - 1 ;NNr 2^(Q1)*r ;% Preallocation of some variablesy zeros(Probl.problem_size,N1) ;fy zeros(Probl.problem_size,N1) ;zn_pred zeros(Probl.problem_size,NNr1) ;if mu 0zn_corr zeros(Probl.problem_size,NNr1) ;elsezn_corr 0 ;end% Evaluation of coefficients of the PECE methodnvett 0 : NNr1 ;nalpha nvett.^alpha ;nalpha1 nalpha.*nvett ;PC.bn nalpha(2:end) - nalpha(1:end-1) ;PC.an [ 1 , (nalpha1(1:end-2) - 2*nalpha1(2:end-1) nalpha1(3:end)) ] ;PC.a0 [ 0 , nalpha1(1:end-2)-nalpha(2:end-1).*(nvett(2:end-1)-alpha-1)] ;PC.halpha1 h^alpha/gamma(alpha1) ;PC.halpha2 h^alpha/gamma(alpha2) ;PC.mu mu ; PC.mu_tol mu_tol ;% Initializing solution and proces of computationt t0 (0 : N)*h ;y(:,1) y0(:,1) ;fy(:,1) f_temp ;[y, fy] Triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl ) ;% Main process of computation by means of the FFT algorithmff [0 2 ] ; nx0 0 ; ny0 0 ;for q 0 : QL 2^q ;[y, fy] DisegnaBlocchi(L, ff, r, Nr, nx0L*r, ny0, t, y, fy, ...zn_pred, zn_corr, N, PC, Probl ) ;ff [ff ff] ; ff(end) 4*L ;end% Evaluation solution in TFINAL when TFINAL is not in the meshif tfinal t(N1)c (tfinal - t(N))/h ;t(N1) tfinal ;y(:,N1) (1-c)*y(:,N) c*y(:,N1) ;endt t(1:N1) ; y y(:,1:N1) ;end% % % r : dimension of the basic square or triangle% L : factor of resizing of the squaresfunction [y, fy] DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, ...zn_pred, zn_corr, N , PC, Probl)nxi nx0 ; nxf nx0 L*r - 1 ;nyi ny0 ; nyf ny0 L*r - 1 ;is 1 ;s_nxi(is) nxi ; s_nxf(is) nxf ; s_nyi(is) nyi ; s_nyf(is) nyf ;i_triangolo 0 ; stop 0 ;while ~stopstop nxir-1 nx0L*r-1 | (nxir-1Nr-1) ; % Ci si ferma quando il triangolo attuale finisce alla fine del quadrato[zn_pred, zn_corr] Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl) ;[y, fy] Triangolo(nxi, nxir-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl) ;i_triangolo i_triangolo 1 ;if ~stopif nxir-1 nxf % Il triangolo finisce dove finisce il quadrato - si scende di livelloi_Delta ff(i_triangolo) ;Delta i_Delta*r ;nxi s_nxf(is)1 ; nxf s_nxf(is) Delta ;nyi s_nxf(is) - Delta 1; nyf s_nxf(is) ;s_nxi(is) nxi ; s_nxf(is) nxf ; s_nyi(is) nyi ; s_nyf(is) nyf ;else % Il triangolo finisce prima del quadrato - si fa un quadrato accantonxi nxi r ; nxf nxi r - 1 ; nyi nyf 1 ; nyf nyf r ;is is 1 ;s_nxi(is) nxi ; s_nxf(is) nxf ; s_nyi(is) nyi ; s_nyf(is) nyf ;endendendend% % function [zn_pred, zn_corr] Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl)coef_beg nxi-nyf ; coef_end nxf-nyi1 ;funz_beg nyi1 ; funz_end nyf1 ;% Evaluation convolution segment for the predictorvett_coef PC.bn(coef_beg:coef_end) ;vett_funz [fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg1) ] ;zzn_pred real(FastConv(vett_coef,vett_funz)) ;zn_pred(:,nxi1:nxf1) zn_pred(:,nxi1:nxf1) zzn_pred(:,nxf-nyf1-1:end-1) ;% Evaluation convolution segment for the correctorif PC.mu 0vett_coef PC.an(coef_beg:coef_end) ;if nyi 0 % Evaluation of the lowest squarevett_funz [zeros(Probl.problem_size,1) , fy(:,funz_beg1:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg1) ] ;else % Evaluation of any square but not the lowestvett_funz [ fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg1) ] ;endzzn_corr real(FastConv(vett_coef,vett_funz)) ;zn_corr(:,nxi1:nxf1) zn_corr(:,nxi1:nxf1) zzn_corr(:,nxf-nyf1:end) ;elsezn_corr 0 ;endend% % function [y, fy] Triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, PC, Probl)for n nxi : min(N,nxf)% Evaluation of the predictorPhi zeros(Probl.problem_size,1) ;if nxi 1 % Case of the first trianglej_beg 0 ;else % Case of any triangle but not the firstj_beg nxi ;endfor j j_beg : n-1Phi Phi PC.bn(n-j)*fy(:,j1) ;endSt StartingTerm(t(n1),Probl.ic) ;y_pred St PC.halpha1*(zn_pred(:,n1) Phi) ;f_pred f_vectorfield(t(n1),y_pred,Probl) ;if PC.mu 0y(:,n1) y_pred ;fy(:,n1) f_pred ;elsej_beg nxi ;Phi zeros(Probl.problem_size,1) ;for j j_beg : n-1Phi Phi PC.an(n-j1)*fy(:,j1) ;endPhi_n St PC.halpha2*(PC.a0(n1)*fy(:,1) zn_corr(:,n1) Phi) ;yn0 y_pred ; fn0 f_pred ;stop 0 ; mu_it 0 ;while ~stopyn1 Phi_n PC.halpha2*fn0 ;mu_it mu_it 1 ;if PC.mu Infstop norm(yn1-yn0,inf) PC.mu_tol ;if mu_it 100 ~stopwarning(MATLAB:fde12:NonConvegence,...strcat(It has been requested to run corrector iterations until convergence but , ...the process does not converge to the tolerance %e in 100 iteration),PC.mu_tol) ;stop 1 ;endelsestop mu_it PC.mu ;endfn1 f_vectorfield(t(n1),yn1,Probl) ;yn0 yn1 ; fn0 fn1 ;endy(:,n1) yn1 ;fy(:,n1) fn1 ;endendend% % function z FastConv(x, y)r length(x) ; problem_size size(y,1) ;z zeros(problem_size,r) ;X fft(x,r) ;for i 1 : problem_sizeY fft(y(i,:),r) ;Z X.*Y ;z(i,:) ifft(Z,r) ;endend% % function f f_vectorfield(t,y,Probl)if isempty(Probl.param)f feval(Probl.fdefun,t,y) ;elsef feval(Probl.fdefun,t,y,Probl.param) ;endend% % function ys StartingTerm(t,ic)ys zeros(size(ic.y0,1),1) ;for k 1 : ic.m_alphays ys (t-ic.t0)^(k-1)/ic.m_alpha_factorial(k)*ic.y0(:,k) ;endend 参考文献更多免费数学建模和仿真教程关注领取
【数据分析】基于数据驱动的分数阶混沌系统建模 附matlab代码
发布时间：2026/5/15 22:45:25
✅作者简介热爱科研的Matlab仿真开发者擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。完整代码获取 定制创新 论文复现点击Matlab科研工作室 关注我领取海量matlab电子书和数学建模资料个人信条做科研博学之、审问之、慎思之、明辨之、笃行之是为博学慎思明辨笃行。 内容介绍一、引言在科学与工程领域混沌系统广泛存在从气象学中的天气预测到电子电路中的信号处理都能发现其身影。传统整数阶混沌系统已得到深入研究但分数阶混沌系统由于其记忆性和遗传性能更精确地描述许多复杂现象。基于数据驱动的分数阶混沌系统建模通过对实际观测数据的分析构建分数阶混沌模型为理解和预测复杂动态系统提供了新的视角和方法。二、分数阶混沌系统基础二分数阶混沌系统特点与整数阶混沌系统相比分数阶混沌系统具有独特的动力学特性。由于分数阶导数的记忆性分数阶混沌系统对初始条件的敏感性更强其相空间轨迹更加复杂。例如在一些分数阶混沌电路中通过调整分数阶的阶数可以观察到不同程度的混沌行为且系统的吸引子形态也会发生变化。这种灵活性为模拟复杂的自然现象和工程系统提供了更强大的工具。三、数据驱动建模的必要性与优势一传统建模方法的局限传统的混沌系统建模通常基于物理原理和数学推导需要对系统的内部结构和参数有深入了解。然而在许多实际场景中系统的物理机制可能非常复杂难以精确描述或者系统参数难以直接测量。例如在生物系统中由于生物过程的高度复杂性和不确定性基于物理原理的建模方法往往无法准确刻画系统的动态行为。二数据驱动建模优势无需精确物理模型数据驱动建模方法直接从观测数据出发无需预先知道系统的物理模型。通过对大量数据的分析和挖掘能够自动发现数据中的潜在规律构建系统模型。这使得数据驱动建模适用于各种复杂系统尤其是那些物理机制不明确的系统。适应复杂动态变化实际系统往往处于动态变化中数据驱动建模可以根据新的数据不断更新模型适应系统的变化。例如在电力系统中负荷需求会随着时间、季节等因素变化基于数据驱动的建模方法能够实时调整模型准确预测电力负荷的变化。四、基于数据驱动的分数阶混沌系统建模步骤一数据采集确定数据源根据研究对象确定合适的数据采集方式和数据源。对于物理实验系统可以通过传感器实时采集系统的输出数据如在混沌电路实验中采集电路的电压、电流等信号。对于一些实际应用场景如气象数据可以从气象站、卫星观测等渠道获取历史数据。保证数据质量采集到的数据应具有足够的长度和精度以反映系统的动态特性。同时要对数据进行预处理去除噪声和异常值。可以采用滤波算法如卡尔曼滤波对数据进行平滑处理提高数据的质量。二特征提取混沌特征分析利用混沌理论中的方法如相空间重构、李雅普诺夫指数计算等分析采集数据的混沌特性。相空间重构通过延迟坐标法将一维时间序列数据映射到高维空间恢复系统的动力学信息。李雅普诺夫指数则用于衡量系统对初始条件的敏感性正的李雅普诺夫指数是混沌系统的重要特征之一。分数阶特征挖掘从数据中挖掘与分数阶相关的特征。由于分数阶系统具有记忆性数据的自相关性在分数阶特征提取中具有重要作用。可以通过分析数据的自相关函数在不同时间尺度上的衰减特性提取与分数阶相关的信息。三模型构建选择建模方法常见的数据驱动建模方法包括神经网络、支持向量机回归等。神经网络具有强大的非线性拟合能力能够逼近任意复杂的函数关系。支持向量机回归则在小样本数据建模中表现出色能够有效避免过拟合问题。根据数据特点和建模需求选择合适的建模方法。确定模型结构以神经网络为例需要确定网络的层数、每层的神经元个数等结构参数。可以通过交叉验证、网格搜索等方法优化模型结构提高模型的泛化能力。对于分数阶混沌系统建模还需要考虑如何将分数阶微积分引入模型中例如通过在神经网络的隐藏层中加入分数阶微分算子构建分数阶神经网络模型。四模型验证与优化模型验证将采集的数据划分为训练集、验证集和测试集。使用训练集对模型进行训练利用验证集调整模型的超参数最后在测试集上验证模型的性能。常用的验证指标包括均方误差MSE、平均绝对误差MAE等这些指标用于衡量模型预测值与实际值之间的偏差。模型优化根据验证结果对模型进行优化。如果模型在测试集上的误差较大可以尝试调整模型结构、增加数据量或采用更复杂的特征提取方法。同时也可以结合其他优化算法如遗传算法、粒子群优化算法等对模型的参数进行全局优化提高模型的性能。⛳️ 运行结果 部分代码function [t, y] fde12(alpha,fdefun,t0,tfinal,y0,h,param,mu,mu_tol)%FDE12 Solves an initial value problem for a non-linear differential% equation of fractional order (FDE). The code implements the% predictor-corrector PECE method of Adams-Bashforth-Moulton type% described in [1].%% [T,Y] FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,h) integrates the initial value% problem for the FDE, or the system of FDEs, of order ALPHA 0% D^ALPHA Y(t) FDEFUN(T,Y(T))% Y^(k)(T0) Y0(:,k1), k0,...,m-1% where m is the smallest integer grater than ALPHA and D^ALPHA is the% fractional derivative according to the Caputos definition. FDEFUN is a% function handle corresponding to the vector field of the FDE and for a% scalar T and a vector Y, FDEFUN(T,Y) must return a column vector. The% set of initial conditions Y0 is a matrix with a number of rows equal to% the size of the problem (hence equal to the number of rows of the% output of FDEFUN) and a number of columns depending on ALPHA and given% by m. The step-size H0 is assumed constant throughout the integration.%% [T,Y] FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM) solves as above with% the additional set of parameters for the FDEFUN as FDEFUN(T,Y,PARAM).%% [T,Y] FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU) solves the FDE with% the selected number MU of multiple corrector iterations. The following% values for MU are admissible:% MU 0 : the corrector is not evaluated and the solution is provided% just by the predictor method (the first order rectangular rule);% MU 0 : the corrector is evaluated by the selected number MU of times;% the classical PECE method is obtained for MU1;% MU Inf : the corrector is evaluated for a certain number of times% until convergence of the iterations is reached (for convergence the% difference between two consecutive iterates is tested).% The defalut value for MU is 1%% [T,Y] FDE12(ALPHA,FDEFUN,T0,TFINAL,Y0,H,PARAM,MU,MU_TOL) allows to% specify the tolerance for testing convergence when MU Inf. If not% specified, the default value MU_TOL 1.E-6 is used.%% FDE12 is an implementation of the predictor-corrector method of% Adams-Bashforth-Moulton studied in [1]. Convergence and accuracy of% the method are studied in [2]. The implementation with multiple% corrector iterations has been proposed and discussed for multiterm FDEs% in [3]. In this implementation the discrete convolutions are evaluated% by means of the FFT algorithm described in [4] allowing to keep the% computational cost proportional to N*log(N)^2 instead of N^2 as in the% classical implementation; N is the number of time-point in which the% solution is evaluated, i.e. N (TFINAL-T)/H. The stability properties% of the method implemented by FDE12 have been studied in [5].%% [1] K. Diethelm, A.D. Freed, The Frac PECE subroutine for the numerical% solution of differential equations of fractional order, in: S. Heinzel,% T. Plesser (Eds.), Forschung und Wissenschaftliches Rechnen 1998,% Gessellschaft fur Wissenschaftliche Datenverarbeitung, Gottingen, 1999,% pp. 57-71.%% [2] K. Diethelm, N.J. Ford, A.D. Freed, Detailed error analysis for a% fractional Adams method, Numer. Algorithms 36 (1) (2004) 31-52.%% [3] K. Diethelm, Efficient solution of multi-term fractional% differential equations using P(EC)mE methods, Computing 71 (2003), pp.% 305-319.%% [4] E. Hairer, C. Lubich, M. Schlichte, Fast numerical solution of% nonlinear Volterra convolution equations, SIAM J. Sci. Statist. Comput.% 6 (3) (1985) 532-541.%% [5] R. Garrappa, On linear stability of predictor-corrector algorithms% for fractional differential equations, Internat. J. Comput. Math. 87% (10) (2010) 2281-2290.%% Copyright (c) 2011-2012, Roberto Garrappa, University of Bari, Italy% garrappa at dm dot uniba dot it% Revision: 1.2 - Date: July, 6 2012% Check inputsif nargin 9mu_tol 1.0e-6 ;if nargin 8mu 1 ;if nargin 7param [] ;endendend% Check order of the FDEif alpha 0error(MATLAB:fde12:NegativeOrder,...[The order ALPHA of the FDE must be positive. The value ...ALPHA %f can not be accepted. See FDE12.], alpha);end% Check the step--size of the methodif h 0error(MATLAB:fde12:NegativeStepSize,...[The step-size H for the method must be positive. The value ...H %e can not be accepted. See FDE12.], h);end% Structure for storing initial conditionsic.t0 t0 ;ic.y0 y0 ;ic.m_alpha ceil(alpha) ;ic.m_alpha_factorial factorial(0:ic.m_alpha-1) ;% Structure for storing information on the problemProbl.ic ic ;Probl.fdefun fdefun ;Probl.problem_size size(y0,1) ;Probl.param param ;% Check number of initial conditionsif size(y0,2) ic.m_alphaerror(MATLAB:fde12:NotEnoughInputs, ...[A not sufficient number of assigned initial conditions. ...Order ALPHA %f requires %d initial conditions. See FDE12.], ...alpha,ic.m_alpha);end% Check compatibility size of the problem with size of the vector fieldf_temp f_vectorfield(t0,y0(:,1),Probl) ;if Probl.problem_size ~ size(f_temp,1)error(MATLAB:fde12:SizeNotCompatible, ...[Size %d of the problem as obtained from initial conditions ...(i.e. the number of rows of Y0) not compatible with the ...size %d of the output of the vector field FDEFUN. ...See FDE12.], Probl.problem_size,size(f_temp,1));end% Number of points in which to evaluate weights and solutionr 16 ;N ceil((tfinal-t0)/h) ;Nr ceil((N1)/r)*r ;Q ceil(log2(Nr/r)) - 1 ;NNr 2^(Q1)*r ;% Preallocation of some variablesy zeros(Probl.problem_size,N1) ;fy zeros(Probl.problem_size,N1) ;zn_pred zeros(Probl.problem_size,NNr1) ;if mu 0zn_corr zeros(Probl.problem_size,NNr1) ;elsezn_corr 0 ;end% Evaluation of coefficients of the PECE methodnvett 0 : NNr1 ;nalpha nvett.^alpha ;nalpha1 nalpha.*nvett ;PC.bn nalpha(2:end) - nalpha(1:end-1) ;PC.an [ 1 , (nalpha1(1:end-2) - 2*nalpha1(2:end-1) nalpha1(3:end)) ] ;PC.a0 [ 0 , nalpha1(1:end-2)-nalpha(2:end-1).*(nvett(2:end-1)-alpha-1)] ;PC.halpha1 h^alpha/gamma(alpha1) ;PC.halpha2 h^alpha/gamma(alpha2) ;PC.mu mu ; PC.mu_tol mu_tol ;% Initializing solution and proces of computationt t0 (0 : N)*h ;y(:,1) y0(:,1) ;fy(:,1) f_temp ;[y, fy] Triangolo(1, r-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl ) ;% Main process of computation by means of the FFT algorithmff [0 2 ] ; nx0 0 ; ny0 0 ;for q 0 : QL 2^q ;[y, fy] DisegnaBlocchi(L, ff, r, Nr, nx0L*r, ny0, t, y, fy, ...zn_pred, zn_corr, N, PC, Probl ) ;ff [ff ff] ; ff(end) 4*L ;end% Evaluation solution in TFINAL when TFINAL is not in the meshif tfinal t(N1)c (tfinal - t(N))/h ;t(N1) tfinal ;y(:,N1) (1-c)*y(:,N) c*y(:,N1) ;endt t(1:N1) ; y y(:,1:N1) ;end% % % r : dimension of the basic square or triangle% L : factor of resizing of the squaresfunction [y, fy] DisegnaBlocchi(L, ff, r, Nr, nx0, ny0, t, y, fy, ...zn_pred, zn_corr, N , PC, Probl)nxi nx0 ; nxf nx0 L*r - 1 ;nyi ny0 ; nyf ny0 L*r - 1 ;is 1 ;s_nxi(is) nxi ; s_nxf(is) nxf ; s_nyi(is) nyi ; s_nyf(is) nyf ;i_triangolo 0 ; stop 0 ;while ~stopstop nxir-1 nx0L*r-1 | (nxir-1Nr-1) ; % Ci si ferma quando il triangolo attuale finisce alla fine del quadrato[zn_pred, zn_corr] Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl) ;[y, fy] Triangolo(nxi, nxir-1, t, y, fy, zn_pred, zn_corr, N, PC, Probl) ;i_triangolo i_triangolo 1 ;if ~stopif nxir-1 nxf % Il triangolo finisce dove finisce il quadrato - si scende di livelloi_Delta ff(i_triangolo) ;Delta i_Delta*r ;nxi s_nxf(is)1 ; nxf s_nxf(is) Delta ;nyi s_nxf(is) - Delta 1; nyf s_nxf(is) ;s_nxi(is) nxi ; s_nxf(is) nxf ; s_nyi(is) nyi ; s_nyf(is) nyf ;else % Il triangolo finisce prima del quadrato - si fa un quadrato accantonxi nxi r ; nxf nxi r - 1 ; nyi nyf 1 ; nyf nyf r ;is is 1 ;s_nxi(is) nxi ; s_nxf(is) nxf ; s_nyi(is) nyi ; s_nyf(is) nyf ;endendendend% % function [zn_pred, zn_corr] Quadrato(nxi, nxf, nyi, nyf, fy, zn_pred, zn_corr, PC, Probl)coef_beg nxi-nyf ; coef_end nxf-nyi1 ;funz_beg nyi1 ; funz_end nyf1 ;% Evaluation convolution segment for the predictorvett_coef PC.bn(coef_beg:coef_end) ;vett_funz [fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg1) ] ;zzn_pred real(FastConv(vett_coef,vett_funz)) ;zn_pred(:,nxi1:nxf1) zn_pred(:,nxi1:nxf1) zzn_pred(:,nxf-nyf1-1:end-1) ;% Evaluation convolution segment for the correctorif PC.mu 0vett_coef PC.an(coef_beg:coef_end) ;if nyi 0 % Evaluation of the lowest squarevett_funz [zeros(Probl.problem_size,1) , fy(:,funz_beg1:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg1) ] ;else % Evaluation of any square but not the lowestvett_funz [ fy(:,funz_beg:funz_end) , zeros(Probl.problem_size,funz_end-funz_beg1) ] ;endzzn_corr real(FastConv(vett_coef,vett_funz)) ;zn_corr(:,nxi1:nxf1) zn_corr(:,nxi1:nxf1) zzn_corr(:,nxf-nyf1:end) ;elsezn_corr 0 ;endend% % function [y, fy] Triangolo(nxi, nxf, t, y, fy, zn_pred, zn_corr, N, PC, Probl)for n nxi : min(N,nxf)% Evaluation of the predictorPhi zeros(Probl.problem_size,1) ;if nxi 1 % Case of the first trianglej_beg 0 ;else % Case of any triangle but not the firstj_beg nxi ;endfor j j_beg : n-1Phi Phi PC.bn(n-j)*fy(:,j1) ;endSt StartingTerm(t(n1),Probl.ic) ;y_pred St PC.halpha1*(zn_pred(:,n1) Phi) ;f_pred f_vectorfield(t(n1),y_pred,Probl) ;if PC.mu 0y(:,n1) y_pred ;fy(:,n1) f_pred ;elsej_beg nxi ;Phi zeros(Probl.problem_size,1) ;for j j_beg : n-1Phi Phi PC.an(n-j1)*fy(:,j1) ;endPhi_n St PC.halpha2*(PC.a0(n1)*fy(:,1) zn_corr(:,n1) Phi) ;yn0 y_pred ; fn0 f_pred ;stop 0 ; mu_it 0 ;while ~stopyn1 Phi_n PC.halpha2*fn0 ;mu_it mu_it 1 ;if PC.mu Infstop norm(yn1-yn0,inf) PC.mu_tol ;if mu_it 100 ~stopwarning(MATLAB:fde12:NonConvegence,...strcat(It has been requested to run corrector iterations until convergence but , ...the process does not converge to the tolerance %e in 100 iteration),PC.mu_tol) ;stop 1 ;endelsestop mu_it PC.mu ;endfn1 f_vectorfield(t(n1),yn1,Probl) ;yn0 yn1 ; fn0 fn1 ;endy(:,n1) yn1 ;fy(:,n1) fn1 ;endendend% % function z FastConv(x, y)r length(x) ; problem_size size(y,1) ;z zeros(problem_size,r) ;X fft(x,r) ;for i 1 : problem_sizeY fft(y(i,:),r) ;Z X.*Y ;z(i,:) ifft(Z,r) ;endend% % function f f_vectorfield(t,y,Probl)if isempty(Probl.param)f feval(Probl.fdefun,t,y) ;elsef feval(Probl.fdefun,t,y,Probl.param) ;endend% % function ys StartingTerm(t,ic)ys zeros(size(ic.y0,1),1) ;for k 1 : ic.m_alphays ys (t-ic.t0)^(k-1)/ic.m_alpha_factorial(k)*ic.y0(:,k) ;endend 参考文献更多免费数学建模和仿真教程关注领取
