matlab解常微分方程常用数值解法1：前向欧拉法和改进的欧拉法

小林up

已于 2023-08-13 12:30:01 修改

阅读量2.1k

点赞数 4

分类专栏： MATLAB 文章标签： matlab 常微分方程数值解欧拉法改进欧拉法

于 2023-08-12 11:53:19 首次发布

本文链接：https://blog.csdn.net/subtitle_/article/details/132245917

版权

MATLAB 专栏收录该内容

45 篇文章 28 订阅

订阅专栏

总结和记录一下matlab求解常微分方程常用的数值解法，本文先从欧拉法和改进的欧拉法讲起。

$\frac{d x}{d t}=f(x, t), \quad x\left(t_{0}\right)=x_{0}$

1. 前向欧拉法

前向欧拉法使用了泰勒展开的第一项线性项逼近。

$x\left(t_{0}+h\right)=x\left(t_{0}\right)+h x^{\prime}\left(t_{0}\right)+\frac{1}{2} x^{\prime \prime}(\xi) h^{2}, \quad t_{0}<\xi<t_{0}+h$

$x_{k+1}=x_{k}+h x'_k+O\left(h^{2}\right)=x_{k}+h f\left(x_{k}, t_{k}\right)+O\left(h^{2}\right)$

2. 改进的欧拉法

在原来前向欧拉法的基础上泰勒展开使用了前面两项：
$x_{k+1}=x_{k}+h x^{\prime}_k+\frac{1}{2} h^{2} x_{k}^{\prime \prime}+O\left(h^{3}\right)$

这里使用：

$x_{k}^{\prime \prime}=\frac{x_{k+1}^{\prime}-x_{k}^{\prime}}{h}$

于是我们有：

$x_{k+1}=x_{k}+\frac{h}{2}\left(x_{k}^{\prime}+x_{k+1}^{\prime}\right)+O\left(h^{3}\right)$

也就是：

$x_{k+1}=x_{k}+\frac{h}{2}\left[f\left(x_{k}, t_{k}\right)+f\left(x_{k+1}, t_{k+1}\right)\right]+O\left(h^{3}\right)$

我们怎么计算 $f(x_{k+1},t_{k+1})$ 呢，因为我们还不知道 $x_{k+1}$ 。

对比前向欧拉法，改进欧拉法的右边不使用 $x_{k+1}$ （我们还不知道），但是我们可以用前向欧拉法计算的 $x_{k}+h f\left(x_{k}, t_{k}\right)$ 来代替 $x_{k+1}$ ，于是我们有
$x_{k+1}=x_{k}+\frac{1}{2}\left(k_{1}+k_{2}\right)+O\left(h^{3}\right), \\\text{其中：}k_{1}=h f\left(x_{k}, t_{k}\right), k_{2}=h f\left(x_{k}+h, t_{k}+k_{1}\right)$

对比一下前向欧拉法：

$x_{k+1}=x_{k}+k_{1}+O\left(h^{2}\right), \quad k_{1}=h f\left(x_{k},t_{k} \right)$

例子

$x^{\prime}=x+t, \quad x(0)=1$

clc
clear

% 测试三个不同的步长
test_times = 3;
% 保存时间、差分时间和步长
h_res=ones(1,test_times);
t_res=cell(1,test_times);%时间
tplot_res=cell(1,test_times);%差分的时间，比时间长度少1
% 保存两种数值方法和解析解的计算结果
x_euler_res=cell(1,test_times);
x_modified_res=cell(1,test_times);
x_exact_res=cell(1,test_times);
% 保存误差
diff1_res=cell(1,test_times);
diff2_res=cell(1,test_times);

for i = 1:test_times
% 设置步长间隔和步长数
	h = 1/10^i; n = 10/h;
	% 设置初始条件
	t=zeros(n+1,1); t(1) = 0;
	x_euler=zeros(n+1,1); x_euler(1) = 1;
	x_modified=zeros(n+1,1); x_modified(1) = 1;
	x_exact=zeros(n+1,1); x_exact(1) = 1;
	% 设置前向欧拉法和改进欧拉法误差存放向量（和精确解比较）
	diff1 = zeros(n,1); diff2 = zeros(n,1); tplot = zeros(n,1);
	% 定义微分方程
	f = inline('xx+tt','tt','xx');
	for k = 1:n
		t(k+1) = t(k) + h;
		% 计算解析解
		x_exact(k+1) = 2*exp(t(k+1)) - t(k+1) - 1;
		% 使用前向欧拉法计算
		k1 = h * f(t(k),x_euler(k));
		x_euler(k+1) = x_euler(k) + k1;
		tplot(k) = t(k+1);
		diff1(k) = x_euler(k+1) - x_exact(k+1);
		diff1(k) = abs(diff1(k) / x_exact(k+1));
		% 使用改进欧拉法计算
		k1 = h * f(t(k),x_modified(k));
		k2 = h * f(t(k+1),x_modified(k)+k1);
		x_modified(k+1) = x_modified(k) + 0.5*(k1+k2);
		diff2(k) = x_modified(k+1) - x_exact(k+1);
		diff2(k) = abs(diff2(k) / x_exact(k+1));
	end
	diff1_res{i}=diff1;
	diff2_res{i}=diff2;
	tplot_res{i}=tplot;
	h_res(i)=h;
	x_euler_res{i}=x_euler;
	x_modified_res{i}=x_modified;
	x_exact_res{i}=x_exact;
	t_res{i}=t;
end

figure
% 不同步长下的解析解、前向欧拉法和改进欧拉法的结果
for i=1:test_times
	subplot(2,2,i)
	plot(t_res{i},x_exact_res{i},'k-',t_res{i},x_euler_res{i},'b--',t_res{i},x_modified_res{i},'r:')
	xlabel('t','Fontsize',18)
	ylabel('x','Fontsize',18)
	legend({'Analytical method','Euler method','modified Euler method'},'Location','best')
	legend boxoff;
	title(['h = ',num2str(h_res(i))]);
end

% 相对误差图
subplot(2,2,4)
for i=1:test_times
	semilogy(tplot_res{i},diff1_res{i},'b-',tplot_res{i},diff2_res{i},'r:')
	hold on
	num1 = 0.2*10/h_res(i); num2 = 0.8*10/h_res(i);
	text(3,diff1_res{i}(num1),['h = ',num2str(h_res(i))],'Fontsize',12,...
	'HorizontalAlignment','right',...
	'VerticalAlignment','bottom')
	text(9,diff2_res{i}(num2),['h = ',num2str(h_res(i))],'Fontsize',12,...
	'HorizontalAlignment','right',...
	'VerticalAlignment','bottom')
end
xlabel('t','Fontsize',18)
ylabel('|relative error|','Fontsize',18)
legend({'Euler method','modified Euler method'},'Location','best')
title('Two Euler method''s relative error')
legend boxoff;