【数据分析】基于有限差分法和乘积积分规则求解分数阶多孔介质方程的Python代码 和matlab代码

​✅作者简介热爱科研的Matlab仿真开发者擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。完整代码获取 定制创新 论文复现点击Matlab科研工作室 关注我领取海量matlab电子书和数学建模资料个人信条做科研博学之、审问之、慎思之、明辨之、笃行之是为博学慎思明辨笃行。 内容介绍一、引言在众多科学与工程领域如石油工程、地下水文学以及土壤科学等多孔介质中流体流动的研究至关重要。分数阶多孔介质方程相较于传统整数阶方程能更精准地描述复杂的非局部和记忆效应。然而由于其分数阶导数的特性求解这类方程颇具挑战。有限差分法结合乘积积分规则为解决这一难题提供了有效的途径它能够将分数阶导数离散化进而实现方程的数值求解。二、分数阶多孔介质方程基础一分数阶导数定义分数阶导数是整数阶导数概念的推广常见的定义包括黎曼 - 刘维尔Riemann - Liouville和卡普托Caputo定义。以卡普托定义为例函数 y(t) 的 α 阶分数阶导数0α1定义为⛳️ 运行结果 部分代码%% Setupb -1;xchebfun(x,[0,1]);ychebfun(y,[-b.^2/41e-5,100]);%% Analytical eigenvaluesdisc sqrt(b.^2 4*y);mu1 -b disc;mu2 - mu1 - 2*b;f mu1.*exp(mu1/2)-mu2.*exp(mu2/2);l roots(f);y chebfun(y,[-100,-b.^2/4-1e-5]);disc sqrt(-b.^2 - 4*y);f -b*sin(disc/2) disc.*cos(disc/2);l max([roots(f);l]);if (b.^24*l0)discl sqrt(b.^2 4*l);u exp((-bdiscl)*x/2) - exp((-b-discl)*x/2);elsediscl sqrt(-b.^2 - 4*l);u exp(-x*b/2).*sin(x*disc(l)/2);endu u/norm(u);%% Eigenvalues with chebfunLchebop((u) diff(u,2) b*diff(u),[0,1]);L.lbcdirichlet;L.rbcneumann;[vv,ll]eigs(L,1);vv vv/norm(vv);%% Compare chebfun and analyticaluu exp(x*b/2).*sin(x*discl/2);uu uu/norm(uu);disp([l,ll])norm(u-vv)%% Eigenvalues for different values of b% Define the range of b valuesb_values 0:10:20;num_b length(b_values);% Create a custom colormap that transitions from blue to redcolors [linspace(0, 1, num_b), zeros(num_b, 1), linspace(1, 0, num_b)];% Initialise figurefigure;hold on;xlabel(Real Part);ylabel(Imaginary Part);title(Eigenvalues for Different Values of b);grid on;axis equal;% Loop through values of b with unique colorsfor k 1:num_bb b_values(k);% Define the operator with current value of bL.op (u) diff(u, 2) b * diff(u);% Compute the first 15 eigenvalueseigenvalues eigs(L, 151) 1e-10i;% Plot eigenvalues with the colour from the colormapplot(eigenvalues, x, Color, colors(k, :), DisplayName, sprintf(b %d, b));end%% Show legend to differentiate b valueslegend show;%% Eigenvalues discretei1;for k 2:12m 2^k1;dx 1.0/2^k;M (diag(2*ones(m,1)) - diag(ones(m-1,1),1) - diag(ones(m-1,1),-1))/(dx*dx) b*(diag(ones(m,1)) - diag(ones(m-1,1),-1))/dx;M(1,2) 0;M(end,end-1)-M(end,end);% [v,l] eigs(inv(M),1);[v,l] sipi(M,eye(size(M,2)),1000,1e-10,-1);disp(l)dof(i) m;lambda(i) l;i i1;end%% test m10m 11;dx 0.1;M (diag(2*ones(m,1)) - diag(ones(m-1,1),1) - diag(ones(m-1,1),-1))/(dx*dx) b*(diag(ones(m,1)) - diag(ones(m-1,1),-1))/dx;M(1,2) 0;M(end,end-1)-M(end,end);% [v,l] eigs(inv(M),1);[v,l] sipi(M,eye(size(M,2)),1000,1e-10,-1);disp(l)%% Functionsfunction [b_k,eigenvalue] sipi(L, M, num_iterations, tolerance, tau)% Compute the largest eigenvalue and eigenvector% using the inverse shifted power iteration method.%% Parameters:% L: Input matrix L (double or sparse).% M: Input matrix M (double or sparse).% num_iterations: Maximum number of iterations (default: 1000).% tolerance: Convergence tolerance (default: 1e-10).% tau: Shift parameter (default: 1.0).%% Returns:% eigenvalue: Dominant eigenvalue.% b_k: Corresponding eigenvector.if nargin 3num_iterations 1000;endif nargin 4tolerance 1e-10;endif nargin 5tau 1.0;end% Start with a random vector of 1sn size(L, 2);b_k ones(n, 1);b_k(1) 0.0; % First component is set to zeroeigenvalue_k 0.0;% Construct the shifted matrixmatrix L - tau * M;disp(matrix)% Print the full matrix% disp(full(matrix));for i 1:num_iterations% Solve the linear systemb_k1 matrix \ b_k;% Compute the Rayleigh quotienteigenvalue_k1 (b_k1 * b_k) / (b_k * b_k);fprintf(Iteration: %d, Eigenvalue: %f\n, i, eigenvalue_k1);% Check for convergenceif abs(eigenvalue_k - eigenvalue_k1) tolerancebreak;end% Update the eigenvector and eigenvalueb_k b_k1 / norm(b_k1);eigenvalue_k eigenvalue_k1;end% Compute the final eigenvalueeigenvalue (1.0 / eigenvalue_k tau);end 参考文献更多免费数学建模和仿真教程关注领取