使用微分矩阵运算符求解 ODE

Question

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

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

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

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.')

然后我在将其发布到此处并收到一个鼓励使用单位矩阵的响应后将其更改为此内容，但我似乎仍然不知所措。

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 解的数值近似值，所以您要寻找代表该 ODE 时间 [=15] 解的数值向量=] 到 x=1。这意味着您的边界条件使得解仅在 0 和 1 之间有效。

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

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

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，我们将设置输出F 中的位置都为 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 现在应该包含从 x=0 到 x=1 到 x=1 的 h = 1/100 步长的 ODE 数值近似值。

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

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

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

希望这对您有所帮助，祝您好运！

使用微分矩阵运算符求解 ODE

USE DIFFERENTIAL MATRIX OPERATOR TO SOLVE ODE

matlab

matrix

numerical-methods

differential-equations

differentiation