在matlab中实现有限差分法

Question

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

EDIT 这是脚本，结果是不正确的

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

用微分方程的右半球。请注意，矩阵中的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'); 

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