思路是通过泰勒展开并将对应阶数的系数抵消来得到差分方程系数。
第一种方法
假设三阶微分的向前差分公式为：
v
i
′
=
a
0
v
i
+
a
1
v
i
+
1
+
a
2
v
i
+
2
+
a
3
v
i
+
3
v'_i=a_0v_{i}+a_1v_{i+1}+a_2v_{i+2}+a_3v_{i+3}
vi′=a0vi+a1vi+1+a2vi+2+a3vi+3
根据泰勒公式：
v
i
+
1
=
v
i
+
h
v
i
′
+
h
2
2
v
i
′
′
+
h
3
6
v
i
(
3
)
+
O
(
h
3
)
v_{i+1}=v_i+hv'_i+\dfrac{h^2}{2}v''_i+\dfrac{h^3}{6}v^{(3)}_i+O(h^3)
vi+1=vi+hvi′+2h2vi′′+6h3vi(3)+O(h3)
v
i
+
2
=
v
i
+
2
h
v
i
′
+
2
h
2
v
i
′
′
+
4
h
3
3
v
i
(
3
)
+
O
(
h
3
)
v_{i+2}=v_i+2hv'_i+2h^2v''_i+\dfrac{4h^3}{3}v^{(3)}_i+O(h^3)
vi+2=vi+2hvi′+2h2vi′′+34h3vi(3)+O(h3)
v
i
+
3
=
v
i
+
3
h
v
i
′
+
9
h
2
2
v
i
′
′
+
9
h
3
2
v
i
(
3
)
+
O
(
h
3
)
v_{i+3}=v_i+3hv'_i+\dfrac{9h^2}{2}v''_i+\dfrac{9h^3}{2}v^{(3)}_i+O(h^3)
vi+3=vi+3hvi′+29h2vi′′+29h3vi(3)+O(h3)
将上式代入向前差分公式并对比系数：
a
0
+
a
1
+
a
2
+
a
3
=
0
a_0+a_1+a_2+a_3=0
a0+a1+a2+a3=0
h
a
1
+
2
h
a
2
+
3
h
a
3
=
1
ha_1+2ha_2+3ha_3=1
ha1+2ha2+3ha3=1
1
2
a
1
+
2
a
2
+
9
2
a
3
=
0
\dfrac{1}{2}a_1+2a_2+\dfrac{9}{2}a_3=0
21a1+2a2+29a3=0
1
6
a
1
+
4
3
a
2
+
9
2
a
3
=
0
\dfrac{1}{6}a_1+\dfrac{4}{3}a_2+\dfrac{9}{2}a_3=0
61a1+34a2+29a3=0
可得
a
0
=
−
11
6
h
,
a
1
=
3
h
,
a
2
=
−
3
2
h
,
a
3
=
1
3
h
a_0=-\dfrac{11}{6h},a_1=\dfrac{3}{h},a_2=-\dfrac{3}{2h},a_3=\dfrac{1}{3h}
a0=−6h11,a1=h3,a2=−2h3,a3=3h1
v i f ′ = − 11 v i + 18 v i + 1 − 9 v i + 2 + 2 v i + 3 6 h v'_{if}= \dfrac{-11v_{i}+18v_{i+1}-9v_{i+2}+2v_{i+3}}{6h} vif′=6h−11vi+18vi+1−9vi+2+2vi+3
将-h代入h可得向后差分公式：
v
i
b
′
=
11
v
i
−
18
v
i
−
1
+
9
v
i
−
2
−
2
v
i
−
3
6
h
v'_{ib}= \dfrac{11v_{i}-18v_{i-1}+9v_{i-2}-2v_{i-3}}{6h}
vib′=6h11vi−18vi−1+9vi−2−2vi−3
v
i
′
=
a
0
v
i
−
3
+
a
1
v
i
−
2
+
a
2
v
i
−
1
+
a
3
v
i
+
a
4
v
i
+
1
+
a
5
v
i
+
2
+
a
6
v
i
+
3
v'_i= a_0v_{i-3}+a_1v_{i-2}+a_2v_{i-1}+a_3v_i+a_4v_{i+1}+a_5v_{i+2}+a_6v_{i+3}
vi′=a0vi−3+a1vi−2+a2vi−1+a3vi+a4vi+1+a5vi+2+a6vi+3
第二种方法
考虑三阶中心差分公式：
v i ′ = a 0 v i − 3 + a 1 v i − 2 + a 2 v i − 1 + a 3 v i + a 4 v i + 1 + a 5 v i + 2 + a 6 v i + 3 v'_i=a_0v_{i-3}+a_1v_{i-2}+a_2v_{i-1}+a_3v_i+a_4v_{i+1}+a_5v_{i+2}+a_6v_{i+3} vi′=a0vi−3+a1vi−2+a2vi−1+a3vi+a4vi+1+a5vi+2+a6vi+3
v
i
+
1
=
v
i
+
h
v
i
′
+
h
2
2
v
i
′
′
+
h
3
6
v
i
(
3
)
+
h
4
4
!
v
i
(
4
)
+
h
5
5
!
v
i
(
5
)
+
h
6
6
!
v
i
(
6
)
+
O
(
h
6
)
v_{i+1}=v_i+hv'_i+\dfrac{h^2}{2}v''_i+\dfrac{h^3}{6}v^{(3)}_i+\dfrac{h^4}{4!}v^{(4)}_i+\dfrac{h^5}{5!}v^{(5)}_i+\dfrac{h^6}{6!}v^{(6)}_i+O(h^6)
vi+1=vi+hvi′+2h2vi′′+6h3vi(3)+4!h4vi(4)+5!h5vi(5)+6!h6vi(6)+O(h6)
v
i
+
2
=
∑
k
=
0
6
(
2
h
)
k
k
!
v
i
k
+
O
(
h
6
)
v_{i+2}=\sum_{k=0}^{6}\dfrac{(2h)^k}{k!}v_i^{k}+O(h^6)
vi+2=∑k=06k!(2h)kvik+O(h6)
v
i
+
3
=
∑
k
=
0
6
(
3
h
)
k
k
!
v
i
k
+
O
(
h
6
)
v_{i+3}=\sum_{k=0}^{6}\dfrac{(3h)^k}{k!}v_i^{k}+O(h^6)
vi+3=∑k=06k!(3h)kvik+O(h6)
v
i
−
1
=
∑
k
=
0
6
(
−
h
)
k
k
!
v
i
k
+
O
(
h
6
)
v_{i-1}=\sum_{k=0}^{6}\dfrac{(-h)^k}{k!}v_i^{k}+O(h^6)
vi−1=∑k=06k!(−h)kvik+O(h6)
v
i
−
2
=
∑
k
=
0
6
(
−
2
h
)
k
k
!
v
i
k
+
O
(
h
6
)
v_{i-2}=\sum_{k=0}^{6}\dfrac{(-2h)^k}{k!}v_i^{k}+O(h^6)
vi−2=∑k=06k!(−2h)kvik+O(h6)
v
i
−
3
=
∑
k
=
0
6
(
−
3
h
)
k
k
!
v
i
k
+
O
(
h
6
)
v_{i-3}=\sum_{k=0}^{6}\dfrac{(-3h)^k}{k!}v_i^{k}+O(h^6)
vi−3=∑k=06k!(−3h)kvik+O(h6)
代入微分公式：
a
0
+
a
1
+
a
2
+
a
3
+
a
4
+
a
5
+
a
6
=
0
a_0+a_1+a_2+a_3+a_4+a_5+a_6=0
a0+a1+a2+a3+a4+a5+a6=0
−
3
a
0
−
2
a
1
−
a
2
+
a
4
+
2
a
5
+
3
a
6
=
1
h
-3a_0-2a_1-a_2+a_4+2a_5+3a_6=\frac{1}{h}
−3a0−2a1−a2+a4+2a5+3a6=h1
9
a
0
+
4
a
1
+
a
2
+
a
4
+
4
a
5
+
9
a
6
=
0
9a_0+4a_1+a_2+a_4+4a_5+9a_6=0
9a0+4a1+a2+a4+4a5+9a6=0
−
27
a
0
−
8
a
1
−
a
2
+
a
4
+
8
a
5
+
27
a
6
=
0
-27a_0-8a_1-a_2+a_4+8a_5+27a_6=0
−27a0−8a1−a2+a4+8a5+27a6=0
81
a
0
+
16
a
1
+
a
2
+
a
4
+
16
a
5
+
81
a
6
=
0
81a_0+16a_1+a_2+a_4+16a_5+81a_6=0
81a0+16a1+a2+a4+16a5+81a6=0
−
243
a
0
−
32
a
1
−
a
2
+
a
4
+
32
a
5
+
243
a
6
=
0
-243a_0-32a_1-a_2+a_4+32a_5+243a_6=0
−243a0−32a1−a2+a4+32a5+243a6=0
729
a
0
+
64
a
1
+
a
2
+
a
4
+
64
a
5
+
729
a
6
=
0
729a_0+64a_1+a_2+a_4+64a_5+729a_6=0
729a0+64a1+a2+a4+64a5+729a6=0
解上述方程可得：
v i ′ = − v i − 3 + 9 v i − 2 − 45 v i − 1 + 45 v i + 1 − 9 v i + 2 + 1 v i + 3 60 h v'_i= \dfrac{-v_{i-3}+9v_{i-2}-45v_{i-1}+45v_{i+1}-9v_{i+2}+1v_{i+3}}{60h} vi′=60h−vi−3+9vi−2−45vi−1+45vi+1−9vi+2+1vi+3
# matlab计算n阶中心差分方程系数
n = 1;
A = zeros(2*n+1);
A(1,:) = ones(1,2*n+1);
for i = 1:2*n
for j = 1:n
A(i+1,n+1-j) = (-j)^i;
A(i+1,n+1+j) = j^i;
end
end
iA = inv(A);
format rat;
W = iA(:,2)
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
