Implement finite difference method in matlab

Question

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.

Why not just calculate values for the new timestep into a separate vector? — Oliver Charlesworth
– Oliver Charlesworth, Commented Jan 12, 2015 at 14:22
How is that loop supposed to work when for the first iteration neither y(i-1) (which is y(0), but you can't index with 0 in MATLAB, it uses 1-based indexing) nor y(i+1) (which is y(2)) is defined? — sobek
– sobek, Commented Jan 12, 2015 at 14:28
Look at the first iteration of your loop, your index i will take the value of 1. Then you try to calculate that value out of y(0) and y(2). But y(0) is not defined in MATLAB, because the first value of the any vector in MATLAB has the index 1. — sobek
– sobek, Commented Jan 12, 2015 at 14:30
Can you specify the differential equation you want to solve using this method? — Nemesis
– Nemesis, Commented Jan 12, 2015 at 14:32

Nemesis · Accepted Answer · 2015-01-12 15:24:41Z

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

You're welcome. I hope the code and the approach is clear so far. Basically solving this kind of diff. eq. is easiest using a system of equations as above.

Collectives™ on Stack Overflow

Implement finite difference method in matlab

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related