MATLAB code for: A singularity-free Crank-Nicolson/Karatay scheme for time-fractional diffusion in spheres with application to soybean hydration
Description
MATLAB code for: A singularity-free Crank-Nicolson/Karatay scheme for time-fractional diffusion in spheres with application to soybean hydration. Corresponding author: Esdras P. de Carvalho (epcarvalho@uem.br). Manuscript submitted to Chemical Engineering Science. This package contains the complete MATLAB implementation used to produce the numerical results, tables and figures of the manuscript. All computation and all figure rendering are done natively in MATLAB (no external language is used anywhere in this pipeline). Version 2 update: figure color palette revised to a colorblind-safe (Okabe–Ito) scheme across Figures 1–9; no change to the numerical methods, fitted parameters, or reported results.
Files
Steps to reproduce
Requirements: MATLAB R2021b with the Optimization Toolbox. No other toolboxes or external software are required. Place all files of this dataset in a single folder (e.g. "Matlab/") and create a sibling folder named "Figures/" at the same level, so the default output path "../Figures" used by the scripts resolves correctly. Add the folder to the MATLAB path (or set it as the current folder). Run "fit_soybean_hydration" to fit the three models (numerical, analytical, classical) at all five temperatures and reproduce Tables 3-5 and Figures 1-7. Optional name-value pairs allow changing the mesh (e.g. fit_soybean_hydration('M',60,'N',2000)). Run "additional_analyses" to reproduce the AICc/BIC comparison, the (D, alpha) identifiability diagnostics, the residual/Durbin-Watson analysis, the mesh-convergence check, the Mittag-Leffler truncation check, and the Arrhenius leave-one-temperature-out cross-validation (Tables 7-9, 11, and the underlying data for Figures 8-9). This script internally calls fit_soybean_hydration.m, which must be on the path. Run "make_fig8_matlab" and "make_fig9_matlab" to generate the final, publication-ready versions of Figures 8 and 9 (Fig_SSE_Contours.jpg and Fig_Residual_Analysis.jpg). Run "numerical_convergence_study" to independently verify the theoretical convergence orders O(dr^2) in space and O(dt^(2-alpha)) in time of the Crank-Nicolson + Karatay scheme. Run "Bootstrap_DAlpha" to reproduce the residual-bootstrap confidence intervals for D and alpha (default B = 200 resamples; expect roughly 10-15 minutes total). A cached result for B = 200 is included as Bootstrap_B200.mat, so this step may be skipped if only the existing bootstrap output needs to be inspected. CI_results.csv contains the corresponding 95% confidence-interval summary already computed. Each script prints its numerical results to the console (matching the values reported in the manuscript's tables) and, where applicable, saves the corresponding figure(s) as .jpg/.eps in the Figures/ folder.