Search code examples
matlabdifferential-equations

Implement finite difference method in matlab

I am trying to implement the finite difference method in matlab. I did some calculations and I got that y(i) is a function of y(i-1) and y(i+1), when I know y(1) and y(n+1). However, I don't know how I can implement this so the values of y are updated the right way. I tried using 2 fors, but it's not going to work that way.

EDIT This is the script and the result isn't right

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

The equation is: y''(x) - 4y'(x) + 3y(x) = x * e ^ (2x), y(0) = 4, y(1) = 6

Thanks.

Solution

  • Consider the following code. The central differential quotient is discretized.

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

    The solution region is specified.

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

    As said in my comment, using this method, all points have to be solved simultaneously. Therefore, above numerical approximation of the equation is transformed in a linear system of euqations.

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

    With the rhs of the differential equation. Note that the known values are calculated by 1 in the matrix and the actual value in the solution vector.

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

y = A\Y;

    Having an equation to approximate the first order derivative (see above) you can verify the solution. (note, ddx2 is an own function)

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

    I hope the solution shown below is correct.

    Solution Verification