1

我们被要求在 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，或者究竟是什么？我一无所知！

4

1 回答 1

1

因为这是ODE 解的数值近似，所以您正在寻找一个数值向量，该向量代表该 ODE 解的时间x=0到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)您在代码中发现的来完成的。

一旦你制定了你的Yand F，你可以使用mldivideor\运算符并通过以下方式求解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定义为ODEh = (b - a)/n中n要求解的总点数。

现在，我们必须包括边界条件。这仅仅意味着我们知道 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现在应该包含您对 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您可以根据边界条件看到这一点。

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

于 2015-06-28T19:29:33.793 回答