使用差分矩阵运算符求解ODE [英] USE DIFFERENTIAL MATRIX OPERATOR TO SOLVE ODE
问题描述
我们被要求在MATLAB上定义自己的微分运算符，然后按照一系列步骤进行操作，然后我们应使用微分运算符来解决边值问题:
We were asked to define our own differential operators on MATLAB, and I did it following a series of steps, and then we should use the differential operators to solve a boundary value problem:
-y''+ 2y'-y = x，y(0)= y(1)= 0
我的代码如下，它用于计算这个(一阶和二阶导数)
my code was as follows, it was used to compute this (first and second derivative)
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的响应后，将其更改为此，但是我似乎仍然无处可寻.
then I changed it to this after posting it here and getting one response that encouraged the usage of Identity Matrix, however I still seem to be getting no where.
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，还是确切地定义什么?我无能为力！
I am not sure how to proceed with this, like should I define y and x, or what exactly? I am clueless!
推荐答案
由于这是ODE解决方案的数值近似值，因此您正在寻找一个代表该数值的矢量.从时间x=0到x=1对此ODE的解.这意味着您的边界条件使之有效，因此解仅在0到1之间有效.
Because this is a numerical approximation to the solution of the ODE, you are seeking to find a numerical vector that is representative of the solution to this ODE from time x=0 to x=1. This means that your boundary conditions make it so that the solution is only valid between 0 and 1.
现在这也是反向问题. 在上一篇文章中我们一起做过，您知道输入向量是什么，并进行矩阵向量乘法运算，就可以对该输入向量进行输出导数运算.现在，为您提供了衍生产品的输出，并且您正在寻找什么是原始输入.现在，这涉及求解线性方程组.
Also this is now the reverse problem. In the previous post we did together, you know what the input vector was, and doing a matrix-vector multiplication produced the output derivative operation on that input vector. Now, you are given the output of the derivative and you are now seeking what the original input was. This now involves solving a linear system of equations.
基本上，您的问题现在是这样:
Essentially, your problem is now this:
YX = F
Y是您从矩阵导数运算符中得出的系数，它是一个n x n矩阵，X是ODE的解，它是一个n x 1向量，而F是您与ODE关联的功能，也是一个n x 1向量.在我们的情况下，该值为x.要找到Y，您已经在代码中完成了很多工作.您只需采用每个矩阵运算符(一阶和二阶导数)，然后将它们与适当的符号和比例相加即可尊重ODE的左侧.顺便说一句，您的一阶导数和二阶导数矩阵是正确的.剩下的就是将-y项添加到组合中，这可以通过-eye(n)完成，正如您在代码中所发现的那样.
Y are the coefficients from the matrix derivative operators that you derived, which is a n x n matrix, X would be the solution to the ODE, which is a n x 1 vector and F would be the function you are associating the ODE with, also a n x 1 vector. In our case, that would be x. To find Y, you've pretty much done that already in your code. You simply take each matrix operator (first and second derivative) and you add them together with the proper signs and scales to respect the left-hand side of the ODE. BTW, your first derivative and second derivative matrices are correct. What's left is adding the -y term to the mix, and that is accomplished by -eye(n) as you have found out in your code.
一旦制定了Y和F，就可以使用 mldivide 或\运算符并求解X并通过以下方式获得此线性系统的解:
Once you formulate your Y and F, you can use the mldivide or \ operator and solve for X and get the solution to this linear system via:
X = Y \ F;
以上从本质上解决了由Y和F形成的线性方程组，并将其存储在X中.
The above essentially solves the linear system of equations formed by Y and F and will be stored in X.
您需要做的第一件事是定义一个从x=0到x=1的点的向量. linspace 可能是最合适的，您可以在其中指定多少个我们要点.现在假设100分:
The first thing you need to do is define a vector of points going from x=0 to x=1. linspace is probably the most suitable where you can specify how many points we want. Let's assume 100 points for now:
x = linspace(0,1,100);
因此，在我们的案例中，h仅为1/100.通常，如果要从起点x = a求解到终点x = b，步长h定义为h = (b - a)/n，其中n是要求解的总点数在ODE中.
Therefore, h in our case is just 1/100. In general, if you want to solve from the starting point x = a up to the end point x = b, the step size h is defined as h = (b - a)/n where n is the total number of points you want to solve for in the ODE.
现在，我们必须包括边界条件.这只是意味着我们知道ODE解决方案的开始和结束.这意味着y(0) = y(1) = 0.因此，我们确保Y的第一行仅将第一列设置为1，而Y的最后一行仅将最后一列设置为1，并将输出位置设置为
Now, we have to include the boundary conditions. This simply means that we know the beginning and ending of the solution of the ODE. This means that y(0) = y(1) = 0. As such, we make sure that the first row of Y has only the first column set to 1 and the last row of Y has only the last column set to 1, and we'll set the output position in F to both be 0. This symbolizes that we already know the solution at these points.
因此，您要解决的最终代码就是:
Therefore, your final code to solve is just:
%// 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的h = 1/100步距ODE的数值近似值.
X should now contain your numerical approximation to the ODE in steps of h = 1/100 starting from x=0 up to x=1.
现在让我们看看它是什么样的:
Now let's see what this looks like:
figure;
plot(x, X);
title('Solution to ODE');
xlabel('x'); ylabel('y');
根据边界条件，您可以看到y(0) = y(1) = 0.
You can see that y(0) = y(1) = 0 as per the boundary conditions.
希望这会有所帮助，祝你好运！
Hope this helps, and good luck!
这篇关于使用差分矩阵运算符求解ODE的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
