求解剛性微分方程用什么數(shù)值解方法?答案就在這里。|何必手算?Matlab求解微分


clc
close all
clear
syms y1(t) y2(t)
equ1 = diff(y1) == -0.01*y1-99.99*y2;
equ2 = diff(y2) == -100*y2;
equ = [equ1, equ2];
cond = [y1(0)==2, y2(0)==1];
S = dsolve(equ, cond)
fun = @(t, y)[-0.01*y(1)-99.99*y(2); -100*y(2)];
tspan = [0 400];
TimerValue = tic;
[t45, y45] = ode45(fun, tspan, [2, 1]);
timer45 = toc(TimerValue);
plot(t45, y45, 'LineWidth', 2);
hold on;
TimerValue = tic;
[t23tb, y23tb] = ode23t(fun, tspan, [2, 1]);
timer23tb = toc(TimerValue);
plot(t23tb, y23tb, 'LineWidth', 2);
TimerValue = tic;
[t15s, y15s] = ode15s(fun, tspan, [2, 1]);
timer15s = toc(TimerValue);
plot(t15s, y15s, 'LineWidth', 2);
legend('ode45', 'ode23tb', 'ode15s');
errorfun = @(t, y)sqrt(sum((exp(-100*t)-y(2)).^2)/length(t))+sqrt(sum(((exp(-100*t)+exp(-0.01*t))-y(1)).^2)/length(t));
error45 = [timer45 errorfun(t45, y45)]
error23tb = [timer23tb errorfun(t23tb, y23tb)]
error15s = [timer15s errorfun(t15s, y15s)]