在matlab中实现有限差分法

Question

我正在尝试在 matlab 中实现有限差分法。我做了一些计算，当我知道y(1) 和y(n+1) 时，我得到y(i) 是y(i-1) 和y(i+1) 的函数。但是，我不知道如何实现这一点，因此y 的值会以正确的方式更新。我尝试使用 2 fors，但它不会那样工作。

编辑 这是脚本，结果不对

n = 10;
m = n+1;
h = 1/m;
x = 0:h:1;
y = zeros(m+1,1);
y(1) = 4;
y(m+1) = 6;
s = y;
for i=2:m
   y(i) = y(i-1)*(-1+(-2)*h)+h*h*x(i)*exp(2*x(i)); 
end

for i=m:-1:2
    y(i) = (y(i) + (y(i+1)*(2*h-1)))/(3*h*h-2);
end

等式是： y''(x) - 4y'(x) + 3y(x) = x * e ^ (2x), y(0) = 4, y(1) = 6

谢谢。

Answer 1

考虑以下代码。中心微商离散化。

% Second order diff. equ.
%          y''              -    4*y'                + 3*y    = x*exp(2*x)
% (y(i+1)-2*y(i)+y(i-1))/h^2-4*(y(i+1)-y(i-1))/(2*h) + 3*y(i) = x(i)*exp(2*x(i));

指定解决方案区域。

x = (0:0.01:1)';   % Solution region
h = min(diff(x));  % distance

正如我在评论中所说，使用这种方法，必须同时解决所有问题。因此，将上述方程的数值逼近转化为线性方程组。

% System of equations
% Matrix of coefficients
A = zeros(length(x));
A(1,1) = 1;     % known solu for first point
A(end,end) = 1; % known solu for last point

% y(i)                                                y''  y
A(2:end-1,2:end-1) = A(2:end-1,2:end-1)+diag(repmat(-2/h^2+3,[length(x)-2 1]));
% y(i-1)                                              y''  -4*y'
A(1:end-1,1:end-1) = A(1:end-1,1:end-1)+diag(repmat(1/h^2+4/(2*h),[length(x)-2 1]),-1);
% y(i+1)                                      y''  -4*y'
A(2:end,2:end) = A(2:end,2:end)+diag(repmat(1/h^2-4/(2*h),[length(x)-2 1]),+1);

与微分方程的rhs。请注意，已知值是由矩阵中的1 和解向量中的实际值计算得出的。

Y = x.*exp(2*x);
Y(1) = 4;   % known solu for first point
Y(end) = 6; % known solu for last point

y = A\Y;

有一个近似一阶导数的方程（见上文），您可以验证解。 （注意，ddx2 是一个自己的函数）

f1 = ddx2(x,y);  % first derivative (own function)
f2 = ddx2(x,f1); % second derivative (own function)

figure;
plot(x,y);
saveas(gcf,'solu1','png');

figure;
plot(x,f2-4*f1+3*y,x,x.*exp(2*x),'ko');
ylim([0 10]);
legend('lhs','rhs','Location','nw');
saveas(gcf,'solu2','png');

我希望下面显示的解决方案是正确的。