使用差分矩阵算求解ODE(USE DIFFERENTIAL MATRIX OPERATOR TO S

我们被要求以定义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:

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

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;

The above essentially solves the linear system of equations formed by Y and F and will be stored in X.

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);

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.

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

You can see that y(0) = y(1) = 0 as per the boundary conditions.

---

Hope this helps, and good luck!

文章来源: USE DIFFERENTIAL MATRIX OPERATOR TO SOLVE ODE