我正在尝试为有限差分法设计算法,但我有点困惑。所讨论的ODE是y'' - 5y'+ 10y = 10x,y(0)= 0且y(1)= 100。所以我需要一种方法以某种方式获得将从关系中乘以“y_i”的系数:
然后将得到的系数存储到矩阵中,矩阵将是我将通过Gauss-Jordan求解的系统的矩阵。问题归结为如何获得这些系数并将它们移动到矩阵。我想过手工计算系数然后只输入矩阵,但是我需要为0.1,0.001和0.001的步长做这个,所以这里真的不是一个可行的选择。
让我们假设ODE的更一般情况
c1 * y''(x) + c2 * y'(x) + c3 * y(x) + c4 * x = 0
与边界条件
y(0) = lo
y(1) = hi
并且你想为x ∈ [0, 1]
解决这个问题,步长为h = 1 / n
(其中n + 1
是样本数)。我们想要解决yi = y(h * i)
。 yi
的范围从i ∈ [0, n]
。为此,我们想要解决一个线性系统
A y = b
每个内部yi
将强加一个线性约束。因此,我们在n - 1
和A
列中有n - 1
行,对应于未知的yi
。
要设置A
和b
,我们可以简单地在我们未知的yi
上滑动一个窗口(我假设从零开始索引)。
A = 0 //the zero matrix
b = 0 //the zero vector
for i from 1 to n - 1
//we are going to create the constraint for yi and store it in row i-1
//coefficient for yi+1
coeff = c1 / h^2 + c2 / h
if i + 1 < n
A(i - 1, i) = coeff
else
b(i - 1) -= coeff * hi //we already know yi+1
//coefficient for yi
coeff = -2 * c1 / h^2 - c2 / h + c3
A(i - 1, i - 1) = coeff
//coefficient for yi-1
coeff = c1 / h^2
if i - 1 > 0
A(i - 1, i - 2) = coeff
else
b(i - 1) -= coeff * lo //we already know yi-1
//coefficient for x
b(i - 1) -= c4 * i * h
next