数学建模社区-数学中国
标题:
用 Euler 法、Heun 法和改进的 Euler 法求解常微分方程
[打印本页]
作者:
2744557306
时间:
2023-12-31 15:58
标题:
用 Euler 法、Heun 法和改进的 Euler 法求解常微分方程
这是一个 MATLAB 函数，用 Euler 法、Heun 法和改进的 Euler 法求解常微分方程，并比较它们的结果。以下是对你的代码的主要部分的解释：
9 r' `. D2 B) K5 `, ]$ x4 S
function y = Euler1(a, b, N, af)
3 O, F0 v. |8 z' h$ ?
h = (b - a) / N;
' S0 B( z! K5 c- {
x(1) = a;
6 A& a. W4 q/ l4 e
y(1) = af;
! W: d* H" @% P1 Y0 u9 \2 @
yg(1) = af;
/ A9 ]9 G+ M' y
yh(1) = af;
0 F- C2 f) Y! @4 I
jqj(1) = af;
: m, S' q7 f1 F5 T
% q# o4 c. P! n6 U. M# W6 y0 L6 v
% 迭代计算
2 q' L- A6 A9 A2 l1 m$ A/ }
for i = 2:N+1
- s) B' B: K, v) {/ X
% Euler 法
$ Q9 i6 H5 u, |+ P4 `
y(i) = y(i-1) + h * f(x(i-1), y(i-1));
" P: h' ?* ^& ^6 v$ j/ A1 j3 f+ Q
# w0 ]1 P% E! p& Z
% Heun 法
. X6 g9 [: {, f* M* O3 I
yh(i) = yh(i-1) + (h/4) * (f(x(i-1), yh(i-1)) + 3 * f(x(i-1) + 2 * h/3, yh(i-1) + 2 * h * f(x(i-1), yh(i-1)) / 3));
2 h) p; B/ d1 X$ Y
/ B# D: m# l# k. q# q
% 改进 Euler 法
8 Q4 g0 s# S7 i& q9 h; W
yg(i) = yg(i-1) + h * (f(x(i-1), y(i-1)) + f(x(i), y(i) + h * f(x(i-1), y(i-1)))) / 2;
2 ~, ?+ U1 g9 d8 H9 n
* m; @& |0 {% w7 _5 \7 C
x(i) = a + (i-1) * h;
- @5 h* _4 W7 O# \& w* K
jqj(i) = x(i) + exp(-x(i));
0 F6 P' R# G8 ~) F* p6 J' z
end
" s; I$ N+ b) n7 }& E
; F: |4 F2 C% S) M
% 计算误差
/ n! v5 I9 }" i
er = sum((y - jqj).^2); % Euler 法误差
. m3 d! l8 z: R1 Q+ L
erg = sum((yg - jqj).^2); % 改进 Euler 法误差
3 h+ _6 i S, u f6 L; p; K
erh = sum((yh - jqj).^2); % Heun 法误差
: i$ S3 @( R8 B* |) ], G0 U
9 L7 T+ M$ B# q* G6 I5 A
% 输出结果和误差
+ X1 u8 O5 W9 b* v" }
[x', y', yg', yh', jqj']
- r, U# y( o" a* z! R2 k
disp(['Euler法误差: ', num2str(er)]);
: R1 q( {+ r, ~
disp(['改进Euler法误差: ', num2str(erg)]);
/ I! N: U9 Y( D" v3 }
disp(['Heun法误差: ', num2str(erh)]);
) e, K* Q- a& c% }
) J# k) r3 ^ y" n% c v
% 绘制结果
8 c4 W. \4 o H7 x) r
plot(x, y, 'r', x, yg, 'b', x, yh, 'k', x, jqj, 'g');
8 w! y8 L1 A' q0 M- Z( t/ t% q1 D; P
legend('Euler法', '改进Euler法', 'Heun法', '精确解');
) j2 F1 `" ?; y9 L- P" w
end
+ _' {3 T0 Z3 N G9 d1 i' T+ F
; ~7 h8 w. |- M6 K# l' G2 g
% 待解的常微分方程的右侧函数
- V1 ^2 J0 f6 L; G: z f x
function result = f(x, y)
0 o, F. {8 i% n
result = -y - exp(-x);
) @4 W* q- i/ S
end
: C: J; x& V1 X% f: f) ?0 D. X
! o! b: ?, z: h' n
这个函数接受四个参数：
% J6 }+ d! F- b7 E ~
" {' Y5 `% M8 [& w4 i
1.a 和 b：求解ODE的时间范围。
$ N* n0 f. x" f* j
2.N：迭代步数。
" ]" R2 ?! d# E- U5 c
3.af：初始条件。
& a) U8 C5 [( |
1 A. h$ | i5 J8 K" ~" r! s
在函数内部，它使用 Euler 法、Heun 法和改进的 Euler 法分别计算解，并计算了它们与精确解的误差。最后，它输出结果矩阵、误差，并通过 plot 函数将结果绘制出来。
8 e+ O7 M& x9 c# R. F- j+ }, ~
你可以通过调用这个函数并提供合适的参数来运行它，例如：
6 t0 \! ?& Q3 k' [; x* E
Euler1(0, 1, 10, 1);
: K/ b5 t8 l+ W5 H$ g9 P
' ` V* P% y. V v( |
这将演示在给定的时间范围内使用三种不同的数值方法求解简单的常微分方程，并比较它们的结果。
4 m# T" Q: \- E/ U, U5 \
' I1 [! s8 B+ l' _* ^
8 b& j) a5 O' m3 D
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5