clear all; clc; close all; clf;
sinan = sinan();
t0 = 0; tfinal = 60; N = 100;
h = (tfinal-t0)/N;
t = t0:h:tfinal;
%% Initial conditions
y0=[100, 10, 20];
y01=[100];
%% Parameters
Lambda = 0.1; beta = 0.4; mu = 0.01; gamma = 0.7; % parameters for model 1
a = 0.2; % parameter for model 2
Models
SIR Epidemic model
subject to
Decay model
subject to ![Initial condition](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAJEAAAAoCAYAAAD+HRieAAAAAXNSR0IArs4c6QAAB9RJREFUeF7tm2fIbdURhh8TMUaxR2JvWAIhapTYsWssxBIjscQWFbFj7x27sSSaEDUQ1GgiaqI/FBRDYkN/aCyIiErsvSAhNmw8l1mX5WGf86299yn3+9gbLlw4a1aZedfMOzPrm4vu6zTQUgNztZTvxDsN0IGoA0FrDXQgaq3CboIORB0GWmtg3CD6EfA58EzrnX9zgm8BmwH/Ad4f8tzddFNoYFwgcp2dgeWB3wWQhm2cpYFjgCuAl4Y9eTdffw2MA0Su8UtgWeCyEQEonfD7wAkdkMYL+XGAaEPgcOBI4K0xHM/1dgNOBP4/hvVGvcRyob/HgJsKFvsOsB2wI/AKsCTMysKvA+4HvhwwRyPZUYNo8QhfVwP/LFDAMIZ8Gzge+Aj4LfDVMCadwByC5xDgQGBR4FfAX6bYh3RBbz8PcGgW1jcArgX+DpzX53I1lh01iPYHfgIcBXw8RkOsCPwhONLTY1x3GEstDBwAmCwsE17IeacCkRf2j8CWwC+Au7PNaGfn9DLroX/TQyvayI602Kgb1YVeBNxTQ7vfBTYF/tUCeHqjc0JRZwFf1Fh/0kPnjpBj2PkBcDNgVjsIRILkCODy8DZe3g96DuLFMhzq4Qz398XvbWRnTTFKT7RDeKA9gdcLLaMC9wPWHoL32iqAtDvw38L157RhqwF/A9aYAkTpwuqFTgPOrQjjXk5D3UER5g35nwZn8rI3kR0piIzJFwASNdPuTwqsI4D2iYMeF665QKzvkJWBvwJnA3e0mWiCsqUg8sKk8CWh7ndeacWlwBORMT8LtJEdKYiWChLoYUT/VJ8x2ZtxLPAG8DPg0amEpvh9oSCTciJD2nQk2KUgSuD4H7A18HAf3RgVbojfHCfNaCM7EEQy9bUirPwY+GGEhj9lm5N3mLZbl7mkpwa0OnBLgGKQF9DF6il0sQtUHNzYrTtv8iVv+L3IVFTwdPtKQDRvEGUzudzDVJ3V8scD8cPRkXxIspvIznYO/TiRxp0P8BC6v3UrCJtGkrzqQfQaGvz52ODGwL+BjYAHCyyXx+vzI64PgwzLD9zLHsA7BftIQxYLT/rTGjJVQ82W2mSmJSDSTlbpzb7qgOj0AF9TWW1fHM6MsdfH+G17QDE/IH/ZBMgJtP8/I1LNJwsMoee7EbCeMSimF0z1jSEnRey3Ym78L/2mE4jyvdYBkQAXCEYXL0td2dmXoyQ7m4r5e9N3ivrDZ2ElQSS4So2XyN1zNYBXAoi6+yiZc5xjSjxR4n7WhuoAwXCk1/996LyurPWmWfYuAZHcR69iaJDn6DY/zDyZiHwRuC3Tbh3juYdT4lb0zt/WYHX20XatUciXgCinAnWAYDhLXFZOWle2VjhTOf08hV7KyQSYWVX66oQzCfVVwF59ahyCzGce9t8MSasE37qmoBjpvgyPOV8bhbFHNWcJiFJh1dBdBwgnRyFY+zWR1YsVcyIH5pwlr5wKFrODPGtzfB1inSsqpZ35/iwWnhk9JMm6nXrjuc1FSf2gdoogMiNxn++NytIjnLcERC5fmqZLL6yd+SXu2Ua2FoiqsqclwgsZinIv5MSlKb5jrWzfDjwUWVT+FiiV6h8JwFhh9RMU9sZ6e0S5PdOeUynChmzpN52ItWfyotwVZZKSYqOPAncFrKG1ka0FohztZmr2aewu286o6ixbm3GcVdRBxcbcFVelw6k4ZiXb0nz6EkhVXCrf9wJkkdiDJYbZrrcQRdMNRHnyY/2nSud5PenPwGHRzW8jWxtEKUTpMSRkVpVtafQ2+pw4gcMnDIPqJAls2wAe3gaimd6rwFNZEc0SQmoYOr8VbksC1kisAVW9ZBRoum651L2F4JnThpWGs7xL7+X1GUhvcTVREpu58k+9v18b2dogSr2oVcMD2SkeVEjcIqrRcqh+DdDkUQyN1qDsnwkY37xYg/INjCGrt2iZPIXy/VokerF9gb0rwu2cBpZ++ykFUbpYevP1err0vUC5soJLpqcgTWSLUvx0wLweIWG1weqj+36f4UTCba8mT//z8YnoyamsV3hbTo0XkAkoK1XUjvK9VFXFBaDKenyaP0zzYt0ZCtPrG6YG9QDNXO0wWO3XG70Wsr6KMAO2s6Dtqv6YobFsSZ0oGT0RVR9KWSt6s+A6S/LMrg7uE/bWj0xrwQCbCkgHHASivCxQBaJ1ounqPpMiC7Y7Rwzx5cPmkT1JfqUEfupFLmPvy2eu7/bZrRdIvfvPy5laWMqa3Q66+I1k64BIzyKadZlupuRzUxdGtzjF4BI5xwzyNoPCmUoTjO4xpbOla3bjGmigFETp9ZtpsjylzrMK3eTF0QaxrVH6pS68xLyUWLtPPZ8vDwyLqSRQumY3roEGSkGkEXWPxtMmf0FhTP558Kg6TzL6pfjp2aiZYp79NV2ngeo6kaSBKhAZDtaMpxMvxLuiJgDo1bLvk3aJtL30r1RXCK4kIczrQe7n1pgvkXYBZOHS8kMdoHZoaKmBKhCZuhuyNPQ/wiA+HCs1/KAtyWX8K4bStz0pPNkxtlTgsxKJprzs7QCWntG6lFncy10Ia4mIBuJVIEoZk9NZ/LMiPUluIei2B36dNWBtg9j26DxOA6MPW6SUEw173W6+GaSBDkQzyJiTOkoHoklpfgat24FoBhlzUkfpQDQpzc+gdTsQzSBjTuooXwOvTitHAleELQAAAABJRU5ErkJggg==)
%% Models
model_1 = @(t, X) [Lambda- beta*X(1)*X(2);
beta*X(1)*X(2)-(gamma+mu)*X(2);
mu*X(2)];
model_2 = @(t, y) -a*y;
% fractional order
alpha = 1;
% Choose a solver
solver_1 = 'caputo';
solver_2 = 'abc';
solver_3 = 'CLASSICAL';
% How to solve a system of equations?
[t_euler_caputo_1, y_euler_caputo_1] = sinan.euler(model_1, alpha, t, h, y0, solver_1);
[t_euler_ABC_1, y_euler_ABC_1] = sinan.euler(model_1, alpha, t, h, y0, solver_2);
[t_euler_clas_1, y_euler_clas_1] = sinan.euler(model_1, alpha, t, h, y0, solver_3);
[t_rk4_caputo_1, y_rk4_caputo_1] = sinan.rk4(model_1, alpha, t, h, y0, solver_1);
[t_rk4_ABC_1, y_rk4_ABC_1] = sinan.rk4(model_1, alpha, t, h, y0, solver_2);
[t_rk4_clas_1, y_rk4_clas_1] = sinan.rk4(model_1, alpha, t, h, y0, solver_3);
% How to solve an equation?
[t_euler_caputo_2, y_euler_caputo_2] = sinan.euler(model_2, alpha, t, h, y01, solver_1);
[t_euler_ABC_2, y_euler_ABC_2] = sinan.euler(model_2, alpha, t, h, y01, solver_2);
[t_euler_clas_2, y_euler_clas_2] = sinan.euler(model_2, alpha, t, h, y01, solver_3);
[t_rk4_caputo_2, y_rk4_caputo_2] = sinan.rk4(model_2, alpha, t, h, y01, solver_1);
[t_rk4_ABC_2, y_rk4_ABC_2] = sinan.rk4(model_2, alpha, t, h, y01, solver_2);
[t_rk4_clas_2, y_rk4_clas_2] = sinan.rk4(model_2, alpha, t, h, y01, solver_3);
figure(1)
plot(t_euler_caputo_1, y_euler_caputo_1(:, 1), 'r', t_euler_ABC_1, y_euler_ABC_1(:, 1), 'b', t_euler_clas_1, y_euler_clas_1(:, 1), 'k')
legend('Euler caputo', 'Euler abc', 'Euler classical')
figure(2)
plot(t_rk4_caputo_1, y_rk4_caputo_1, 'r', t_rk4_ABC_1, y_rk4_ABC_1, 'b', t_rk4_clas_1, y_rk4_clas_1, 'k')
legend('rk4 caputo', 'rk4 abc', 'rk4 classical')
figure(3)
plot(t_euler_clas_1, y_euler_clas_1, 'r', t_rk4_clas_1, y_rk4_clas_1, 'b')
legend('euler classical', 'rk4 classical')
figure(4)
plot(t_euler_caputo_2, y_euler_caputo_2, 'r', t_euler_ABC_2, y_euler_ABC_2, 'b', t_euler_clas_2, y_euler_clas_2, 'k')
legend('Euler caputo', 'Euler abc', 'Euler classical')
xlabel('time (t)')
ylabel('y(t)')
figure(5)
plot(t_rk4_caputo_2, y_rk4_caputo_2, 'r', t_rk4_ABC_2, y_rk4_ABC_2, 'b', t_rk4_clas_2, y_rk4_clas_2, 'k')
legend('rk4 caputo', 'rk4 abc', 'rk4 classical')
xlabel('time (t)')
ylabel('y(t)')
figure(6)
plot(t_euler_clas_2, y_euler_clas_2, 'r', t_rk4_clas_2, y_rk4_clas_2, 'b')
legend('euler classical', 'rk4 classical')
xlabel('time (t)')
ylabel('y(t)')
引用
Muhammad Sinan (2024). sinan (https://www.mathworks.com/matlabcentral/fileexchange/132812-sinan), MATLAB Central File Exchange. 取得済み .
Sinan, Muhammad. Sinan. Zenodo, 2023, doi:10.5281/ZENODO.8189815.
MATLAB リリースの互換性
作成:
R2022b
すべてのリリースと互換性あり
プラットフォームの互換性
Windows macOS Linuxタグ
謝辞
ヒントを得たファイル: Digital Fractional Order Differentiator/integrator - IIR type, 微分方程式論(Introduction of Differential Equations), Non-linear equations system solver (Newton Raphson), Tutorial on solving DDEs with DDE23, Fractional Differential Equation Exact Solution, Tutorial on solving BVPs with BVP4C, Ordinary Differential Equation Toolbox: ODEbox Version 1.1, Qualitative Analysis of ODEs
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!