数学建模社区-数学中国
标题:
四阶RK法求解常微分方程
[打印本页]
作者:
2744557306
时间:
2023-12-31 17:32
标题:
四阶RK法求解常微分方程
这段MATLAB代码实现了四阶Runge-Kutta（RK）方法来求解常微分方程（ODE）。以下是代码的主要解释：
' y7 b& {4 N3 g' @4 M, G4 G; D
function y = RK(a, b, N, af)
0 Y. X& D$ q( U! c0 n B" x
h = (b - a) / N;
8 R1 w0 U" j2 w+ i5 @
x(1) = a;
1 A1 p# g0 v2 }
y(1) = af;
8 W- ]$ x; r/ `% \. v3 T
jqj(1) = af;
) ]' S0 p2 G6 X0 p% W6 ]
1 `6 G5 B, E& O5 b1 n4 @1 E0 \. v
for i = 2:N+1
$ q: \( U" F* D( |. I
K1 = f(x(i-1), y(i-1));
9 J9 e; g E, @
K2 = f(x(i-1) + h/2, y(i-1) + h*K1/2);
2 b8 ]- q/ s( }# ~4 c7 @4 r9 e& `
K3 = f(x(i-1) + h/2, y(i-1) + h*K2/2);
9 }) z" v5 z8 i, [! k7 f, z9 W8 g
K4 = f(x(i-1) + h, y(i-1) + h*K3);
! s, S3 a8 o$ A$ l/ r
/ W! h3 [" W, O: E, J
y(i) = y(i-1) + (K1 + 2*K2 + 2*K3 + K4) / 6;
& u3 S0 X4 Y2 x5 N, N/ L
x(i) = x(i-1) + (i-1) * h;
. ^$ w+ k" ^' ~3 u. W. r: g
jqj(i) = x(i) + exp(-x(i));
# ?$ D# i3 @7 ~ e! A$ h0 H- J
end
1 ~8 @! A% U2 d# z7 i- s7 D. I
, X5 W4 ?8 I$ W- t
[x', y', jqj']
/ g! e( U+ J: t
er = norm(y - jqj, 2) / norm(y);
. c* A1 h( p% `. A+ n. [
$ O6 Z$ i' k6 d, r2 X' W! h v
plot(x', y', 'r', x', jqj', 'g');
# v) d$ W0 b$ C; L
legend('RK法', '精确解');
- J4 C, V* g3 }# y6 ]
end
$ D" ]3 [2 A$ q9 f/ @
$ |4 S; G3 K0 T6 N
这是代码的简要说明：
& z5 P) `2 B+ k6 E) p
! U* A$ X7 m/ L5 R/ h7 t
1.函数RK接受初始值和终止值a和b，步数N以及解的初始值af。
+ M3 n" A, ~" o& ~4 m
2.初始化数组x、y和jqj，用于存储自变量、使用RK方法得到的解以及精确解的值。
: c6 u5 D) f7 V6 @, V) S) W# q# M
3.for循环执行RK方法的迭代，每一步更新解y。
! @, o" ^0 Z# P* M$ W' ^
4.与RK解同时计算精确解jqj。
1 ^/ G# r. m. I* D2 W
5.该函数以表格形式打印x、y和jqj的值，计算相对误差er，并绘制RK解和精确解的图表。
2 q6 Y, ^/ ~# m- I: x
6.注意：函数f被假定在您的代码中其他地方已经定义，并且表示要使用RK方法求解的函数的导数。
& M+ ]. i+ B% ^0 a$ c* Y) C
, y) @, D/ `# V6 m, R; e
/ l7 S8 T, M4 I1 t, z" t
5 V& m k: b+ S
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5