数学建模社区-数学中国
标题:
四阶RK法求解常微分方程
[打印本页]
作者:
2744557306
时间:
2023-12-31 17:32
标题:
四阶RK法求解常微分方程
这段MATLAB代码实现了四阶Runge-Kutta（RK）方法来求解常微分方程（ODE）。以下是代码的主要解释：
0 I4 j# c$ g3 z. r
function y = RK(a, b, N, af)
$ z4 k- A2 }: n) m( Y
h = (b - a) / N;
& b' r X0 ^: r
x(1) = a;
- i+ P( ]' G3 H# E. \+ F
y(1) = af;
: ^/ H( r: X9 _& k; L, l$ T5 f6 d6 {
jqj(1) = af;
4 r/ z" O) d7 a; y; r
* U, y, j/ u* w( R$ c$ P* A
for i = 2:N+1
2 E u7 z6 ^' M+ \4 T
K1 = f(x(i-1), y(i-1));
; w- Y$ J2 l- \% w; x
K2 = f(x(i-1) + h/2, y(i-1) + h*K1/2);
- B' l! g* M0 u- |2 C1 V" c
K3 = f(x(i-1) + h/2, y(i-1) + h*K2/2);
5 g- F3 m9 k4 Q: e O1 `' u6 M( ^& e
K4 = f(x(i-1) + h, y(i-1) + h*K3);
$ W) f, w$ q, `2 T$ P. q
- z( K e% }! P( x) j
y(i) = y(i-1) + (K1 + 2*K2 + 2*K3 + K4) / 6;
* U3 G4 v, W q, z% a" |; J
x(i) = x(i-1) + (i-1) * h;
& N- j$ I9 ~; m3 F4 v6 r
jqj(i) = x(i) + exp(-x(i));
1 N1 y. I4 c# e8 N( W% l4 l" ]
end
A6 Y9 u, D* D/ m+ p, c, t& n
4 a/ ~8 |8 K# N/ ?( M' m
[x', y', jqj']
k0 |! |) z6 x! C( D2 d
er = norm(y - jqj, 2) / norm(y);
) m* R6 [4 B; A& `, v, E
5 A7 _; p& {0 K; F3 l( q
plot(x', y', 'r', x', jqj', 'g');
: ?6 ?. w0 O7 r! Q5 j/ o
legend('RK法', '精确解');
2 v. u, c; e+ H" U; z& ^9 m! K5 m5 U
end
3 E, K2 D* `9 F5 Z4 b( r, T. b4 C2 q
/ [% m, {" l. w L1 q
这是代码的简要说明：
( k% D) X5 ~# m6 Z! x9 n* ?7 H/ b7 C
' M& }) ?5 N% n1 |1 M
1.函数RK接受初始值和终止值a和b，步数N以及解的初始值af。
) R) R/ @; x3 ]7 h+ `, Q( |! `
2.初始化数组x、y和jqj，用于存储自变量、使用RK方法得到的解以及精确解的值。
( K/ P( v* O' g1 Y \
3.for循环执行RK方法的迭代，每一步更新解y。
# {8 _/ ^3 [& i% s) Q% K
4.与RK解同时计算精确解jqj。
+ A, X+ `; ? ~
5.该函数以表格形式打印x、y和jqj的值，计算相对误差er，并绘制RK解和精确解的图表。
' V& J2 a0 @- Y+ V! v7 J
6.注意：函数f被假定在您的代码中其他地方已经定义，并且表示要使用RK方法求解的函数的导数。
( ]8 D9 U* ^% }& @. D6 R- e
- R* T7 {' a' ?
' N1 u5 [- e* x8 `, w! n( D
$ ?* _5 ~8 \3 ~( O: ^
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5