使用微分矩阵算子求解 ODE

Question

我们被要求在 MATLAB 上定义自己的微分算子，我按照一系列步骤完成了，然后我们应该使用微分算子来解决边值问题：

-y'' + 2y' - y = x, y(0) = y(1) =0

我的代码如下，用于计算this（一阶和二阶导数）

h = 2;
x = 2:h:50; 
y = x.^2 ;

n=length(x);
uppershift = 1; 
U = diag(ones(n-abs(uppershift),1),uppershift);
lowershift = -1;
L= diag(ones(n-abs(lowershift),1),lowershift);
% the code above creates the upper and lower shift matrix


D  = ((U-L))/(2*h); %first differential operator 
D2 = (full (gallery('tridiag',n)))/ -(h^2); %second differential operator

d1= D*y.'
d2= ((D2)*y.')

然后我在将其发布到此处并收到一个鼓励使用 Identity Matrix 的回复后将其更改为此，但我似乎仍然无处可去。

h = 2;

n=10;
uppershift = 1; 
U = diag(ones(n-abs(uppershift),1),uppershift);
lowershift = -1;
L= diag(ones(n-abs(lowershift),1),lowershift);


D  = ((U-L))/(2*h); %first differential operator 
D2 = (full (gallery('tridiag',n)))/ -(h^2); %second differential operator


I= eye(n);

eqn=(-D2 + 2*D - I)*y == x

solve(eqn,y)

我不确定如何进行此操作，例如我应该定义 y 和 x，或者究竟是什么？我一无所知！

Answer 1

因为这是 ODE 解的数值近似，您正在寻找一个数值向量，该向量表示从时间 x=0 到 x=1 的此 ODE 解。这意味着您的边界条件使解决方案仅在 0 和 1 之间有效。

这也是现在的reverse问题。 In the previous post we did together，您知道输入向量是什么，并且进行矩阵向量乘法会在该输入向量上产生输出导数运算。现在，您得到了导数的输出，您现在正在寻找原始输入。这现在涉及求解线性方程组。

基本上，您现在的问题是：

YX = F

Y 是您导出的矩阵导数运算符的系数，这是一个 n x n 矩阵，X 是 ODE 的解，这是一个 n x 1 向量，F 是与 ODE 关联的函数，也是一个 n x 1 向量。在我们的例子中，这将是x。要找到Y，您已经在代码中完成了很多工作。您只需获取每个矩阵运算符（一阶和二阶导数），然后将它们与适当的符号和比例相加，以尊重 ODE 的左侧。顺便说一句，您的一阶导数和二阶导数矩阵是正确的。剩下的就是添加-y 术语，这是由-eye(n) 完成的，正如您在代码中发现的那样。

一旦你制定了Y 和F，你就可以使用mldivide 或\ 运算符求解X 并通过以下方式得到这个线性系统的解：

X = Y \ F;

上面基本上解决了由Y和F组成的线性方程组，并将存储在X中。

您需要做的第一件事是定义一个从x=0 到x=1 的点向量。 linspace 可能是最合适的，您可以指定我们想要的点数。我们现在假设 100 分：

x = linspace(0,1,100);

因此，在我们的例子中，h 就是1/100。一般来说，如果要从起点x = a求解到终点x = b，则步长h定义为h = (b - a)/n，其中n是要求解的总点数在 ODE 中。

现在，我们必须包括边界条件。这仅仅意味着我们知道 ODE 解的开始和结束。这意味着y(0) = y(1) = 0。因此，我们确保Y 的第一行只有第一列设置为1，Y 的最后一行只有最后一列设置为1，我们将在@987654361 中设置输出位置@ 都为 0。这表示我们已经知道这些点的解。

因此，您最终要解决的代码就是：

%// Setup
a = 0; b = 1; n = 100;
x = linspace(a,b,n);
h = (b-a)/n;

%// Your code
uppershift = 1;
U = diag(ones(n-abs(uppershift),1),uppershift);
lowershift = -1;
L= diag(ones(n-abs(lowershift),1),lowershift);
D  = ((U-L))/(2*h); %first differential operator
D2 = (full (gallery('tridiag',n)))/ -(h^2);

%// New code - Create differential equation matrix
Y = (-D2 + 2*D - eye(n));

%// Set boundary conditions on system
Y(1,:) = 0; Y(1,1) = 1;
Y(end,:) = 0; Y(end,end) = 1;

%// New code - Create F vector and set boundary conditions
F = x.';
F(1) = 0; F(end) = 0;

%// Solve system
X = Y \ F;

X 现在应该包含您对 ODE 的数值近似，步长为 h = 1/100，从 x=0 到 x=1。

现在让我们看看它是什么样子的：

figure; 
plot(x, X);
title('Solution to ODE');
xlabel('x'); ylabel('y');

您可以根据边界条件看到y(0) = y(1) = 0。

---

希望这会有所帮助，祝你好运！