diy网站开发公司,大型彩灯制作公司,星河网站建设,在山东省建设监理协会网站总结和记录一下matlab求解常微分方程常用的数值解法#xff0c;本文先从欧拉法和改进的欧拉法讲起。 d x d t f ( x , t ) , x ( t 0 ) x 0 \frac{d x}{d t}f(x, t), \quad x\left(t_{0}\right)x_{0} dtdxf(x,t),x(t0)x0
1. 前向欧拉法
前向欧拉法使用了泰勒展开的第…总结和记录一下matlab求解常微分方程常用的数值解法本文先从欧拉法和改进的欧拉法讲起。 d x d t f ( x , t ) , x ( t 0 ) x 0 \frac{d x}{d t}f(x, t), \quad x\left(t_{0}\right)x_{0} dtdxf(x,t),x(t0)x0
1. 前向欧拉法
前向欧拉法使用了泰勒展开的第一项线性项逼近。 x ( t 0 h ) x ( t 0 ) h x ′ ( t 0 ) 1 2 x ′ ′ ( ξ ) h 2 , t 0 ξ t 0 h 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}\xit_{0}h x(t0h)x(t0)hx′(t0)21x′′(ξ)h2,t0ξt0h x k 1 x k h x k ′ O ( h 2 ) x k h f ( x k , t k ) O ( h 2 ) x_{k1}x_{k}h x_kO\left(h^{2}\right)x_{k}h f\left(x_{k}, t_{k}\right)O\left(h^{2}\right) xk1xkhxk′O(h2)xkhf(xk,tk)O(h2)
2. 改进的欧拉法
在原来前向欧拉法的基础上泰勒展开使用了前面两项 x k 1 x k h x k ′ 1 2 h 2 x k ′ ′ O ( h 3 ) x_{k1}x_{k}h x^{\prime}_k\frac{1}{2} h^{2} x_{k}^{\prime \prime}O\left(h^{3}\right) xk1xkhxk′21h2xk′′O(h3)
这里使用 x k ′ ′ x k 1 ′ − x k ′ h x_{k}^{\prime \prime}\frac{x_{k1}^{\prime}-x_{k}^{\prime}}{h} xk′′hxk1′−xk′
于是我们有 x k 1 x k h 2 ( x k ′ x k 1 ′ ) O ( h 3 ) x_{k1}x_{k}\frac{h}{2}\left(x_{k}^{\prime}x_{k1}^{\prime}\right)O\left(h^{3}\right) xk1xk2h(xk′xk1′)O(h3)
也就是 x k 1 x k h 2 [ f ( x k , t k ) f ( x k 1 , t k 1 ) ] O ( h 3 ) x_{k1}x_{k}\frac{h}{2}\left[f\left(x_{k}, t_{k}\right)f\left(x_{k1}, t_{k1}\right)\right]O\left(h^{3}\right) xk1xk2h[f(xk,tk)f(xk1,tk1)]O(h3)
我们怎么计算 f ( x k 1 , t k 1 ) f(x_{k1},t_{k1}) f(xk1,tk1)呢因为我们还不知道 x k 1 x_{k1} xk1。
对比前向欧拉法改进欧拉法的右边不使用 x k 1 x_{k1} xk1我们还不知道但是我们可以用前向欧拉法计算的 x k h f ( x k , t k ) x_{k}h f\left(x_{k}, t_{k}\right) xkhf(xk,tk)来代替 x k 1 x_{k1} xk1于是我们有 x k 1 x k 1 2 ( k 1 k 2 ) O ( h 3 ) , 其中 k 1 h f ( x k , t k ) , k 2 h f ( x k h , t k k 1 ) x_{k1}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) xk1xk21(k1k2)O(h3),其中k1hf(xk,tk),k2hf(xkh,tkk1)
对比一下前向欧拉法 x k 1 x k k 1 O ( h 2 ) , k 1 h f ( x k , t k ) x_{k1}x_{k}k_{1}O\left(h^{2}\right), \quad k_{1}h f\left(x_{k},t_{k} \right) xk1xkk1O(h2),k1hf(xk,tk)
例子 x ′ x t , x ( 0 ) 1 x^{\prime}xt, \quad x(0)1 x′xt,x(0)1
clear% 测试三个不同的步长
test_times 3;
% 保存时间、差分时间和步长
h_resones(1,test_times);
t_rescell(1,test_times);%时间
tplot_rescell(1,test_times);%差分的时间比时间长度少1
% 保存两种数值方法和解析解的计算结果
x_euler_rescell(1,test_times);
x_modified_rescell(1,test_times);
x_exact_rescell(1,test_times);
% 保存误差
diff1_rescell(1,test_times);
diff2_rescell(1,test_times);for i 1:test_times
% 设置步长间隔和步长数
h 1/10^i; n 10/h;
% set up initial conditions
tzeros(n1,1); t(1) 0;
x_eulerzeros(n1,1); x_euler(1) 1;
x_modifiedzeros(n1,1); x_modified(1) 1;
x_exactzeros(n1,1); x_exact(1) 1;
% 设置不同的比较误差的图
diff1 zeros(n,1); diff2 zeros(n,1); tplot zeros(n,1);
% define right side of differential equation, Equation 1.7.10
f inline(xxtt,tt,xx);
for k 1:n
t(k1) t(k) h;
% 计算解析解
x_exact(k1) 2*exp(t(k1)) - t(k1) - 1;
% 使用前向欧拉法计算
k1 h * f(t(k),x_euler(k));
x_euler(k1) x_euler(k) k1;
tplot(k) t(k1);
diff1(k) x_euler(k1) - x_exact(k1);
diff1(k) abs(diff1(k) / x_exact(k1));
% 使用改进欧拉法计算
k1 h * f(t(k),x_modified(k));
k2 h * f(t(k1),x_modified(k)k1);
x_modified(k1) x_modified(k) 0.5*(k1k2);
diff2(k) x_modified(k1) - x_exact(k1);
diff2(k) abs(diff2(k) / x_exact(k1));
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;
endfigure
for i1: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(TIME,Fontsize,18)
ylabel(|RELATIVE ERROR|,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 i1: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(TIME,Fontsize,18)
ylabel(|RELATIVE ERROR|,Fontsize,18)
legend({Euler method,modified Euler method},Location,best)
legend boxoff;
我们对各个不同的步长进行了比较并比较了它们的误差发现改进的欧拉法要比前向欧拉法的精度更高。随着步长的变小误差也在变小。
