Matlab 常微分方程数值解法：显式 Euler 法、隐式 Euler 预估-校正法、四阶经典 Runge-Kutta 法

% 显式 Euler 法 clear all; clc; h = 0.05; t = 0:h:1; u(1) = 1; v(1) = 2; for i = 1:length(t)-1 u(i+1) = u(i) + h*(-4u(i)+6v(i)); v(i+1) = v(i) + h*(3u(i)-7v(i)); end exact_u = exp(-4t) + 2exp(-t); exact_v = 2exp(-t) - exp(-7t); error_u = abs(u-exact_u); error_v = abs(v-exact_v); subplot(2,1,1); plot(t,error_u); title('Error of u(t) using Explicit Euler Method'); xlabel('t'); ylabel('error'); subplot(2,1,2); plot(t,error_v); title('Error of v(t) using Explicit Euler Method'); xlabel('t'); ylabel('error');

% 隐式 Euler 预估-校正法 clear all; clc; h = 0.1; t = 0:h:1; u(1) = 1; v(1) = 2; for i = 1:length(t)-1 u_p = u(i) + h*(-4u(i)+6v(i)); % 预估 v_p = v(i) + h*(3u(i)-7v(i)); u(i+1) = u(i) + h*(-4u_p+6v_p); % 校正 v(i+1) = v(i) + h*(3u_p-7v_p); end exact_u = exp(-4t) + 2exp(-t); exact_v = 2exp(-t) - exp(-7t); error_u = abs(u-exact_u); error_v = abs(v-exact_v); subplot(2,1,1); plot(t,error_u); title('Error of u(t) using Implicit Euler Method'); xlabel('t'); ylabel('error'); subplot(2,1,2); plot(t,error_v); title('Error of v(t) using Implicit Euler Method'); xlabel('t'); ylabel('error');

% 四阶经典 Runge-Kutta 法 clear all; clc; h = 0.25; t = 0:h:1; u(1) = 1; v(1) = 2; for i = 1:length(t)-1 k1_u = -4u(i) + 6v(i); k1_v = 3u(i) - 7v(i); k2_u = -4*(u(i)+h/2k1_u) + 6(v(i)+h/2k1_v); k2_v = 3(u(i)+h/2k1_u) - 7(v(i)+h/2k1_v); k3_u = -4(u(i)+h/2k2_u) + 6(v(i)+h/2k2_v); k3_v = 3(u(i)+h/2k2_u) - 7(v(i)+h/2k2_v); k4_u = -4(u(i)+hk3_u) + 6(v(i)+hk3_v); k4_v = 3(u(i)+hk3_u) - 7(v(i)+hk3_v); u(i+1) = u(i) + h/6(k1_u+2k2_u+2k3_u+k4_u); v(i+1) = v(i) + h/6*(k1_v+2k2_v+2k3_v+k4_v); end exact_u = exp(-4t) + 2exp(-t); exact_v = 2exp(-t) - exp(-7t); error_u = abs(u-exact_u); error_v = abs(v-exact_v); subplot(2,1,1); plot(t,error_u); title('Error of u(t) using Classic Runge-Kutta Method'); xlabel('t'); ylabel('error'); subplot(2,1,2); plot(t,error_v); title('Error of v(t) using Classic Runge-Kutta Method'); xlabel('t'); ylabel('error');

% 不同步长误差比较 clear all; clc; h = [0.5, 0.25, 0.125]; for j = 1:length(h) t = 0:h(j):1; u(1) = 1; v(1) = 2; for i = 1:length(t)-1 k1_u = -4u(i) + 6v(i); k1_v = 3u(i) - 7v(i); k2_u = -4*(u(i)+h(j)/2k1_u) + 6(v(i)+h(j)/2k1_v); k2_v = 3(u(i)+h(j)/2k1_u) - 7(v(i)+h(j)/2k1_v); k3_u = -4(u(i)+h(j)/2k2_u) + 6(v(i)+h(j)/2k2_v); k3_v = 3(u(i)+h(j)/2k2_u) - 7(v(i)+h(j)/2k2_v); k4_u = -4(u(i)+h(j)k3_u) + 6(v(i)+h(j)k3_v); k4_v = 3(u(i)+h(j)k3_u) - 7(v(i)+h(j)k3_v); u(i+1) = u(i) + h(j)/6(k1_u+2k2_u+2k3_u+k4_u); v(i+1) = v(i) + h(j)/6*(k1_v+2k2_v+2k3_v+k4_v); end exact_u = exp(-4t) + 2exp(-t); exact_v = 2exp(-t) - exp(-7t); error(j) = norm([u(end)-exact_u(end), v(end)-exact_v(end)]); end log_h = log(h); log_e = log(error); plot(log_h, log_e); title('Overall Error using Classic Runge-Kutta Method'); xlabel('log(h)'); ylabel('log(enum)');