Finite Difference Coding Mistake

1 view (last 30 days)
I cannot understand the reason why the approximation obtained through finite difference does not converge to the exact solution of the following problem
clear;
clc;
f=@(x)2*exp(x).*(2*sin(pi*x)-2*pi*cos(pi*x)+pi^2*sin(pi*x));
sol=@(x)2*exp(x).*sin(pi*x);%exact solution
a=0;
b=1;
N=10000;
h=(b-a)/(N+1);
x=(a:h:b)';
x(1)=[];
x(end)=[];
%define diffusion matrix
Ad=2*diag(ones(N,1));
Ad=Ad-diag(ones(N-1,1),-1);
Ad=Ad-diag(ones(N-1,1),1);
Ad=Ad/h^2;
sigma=3;
%define transport matrix
At=zeros(N,N);
At=At+diag(ones(N-1,1),1);
At=At-diag(ones(N-1,1),-1);
At=At*sigma/(2*h);
A=Ad+At;
F=f(x);
U=A\F;
U=[0 ; U ; 0];
x=[a ; x ; b];
plot(x,U,'--');%plot approximation
hold on;
plot(x,sol(x));%plot exact solution

Accepted Answer

infinity
infinity on 20 Jun 2019
In your code, there was a mistake,
%define transport matrix
% At=zeros(N,N);
% At=At+diag(ones(N-1,1),1);
% At=At-diag(ones(N-1,1),-1);
% At=At*sigma/(2*h);
since you did wrong approximtion of 3u(x). You can look at the code below, it will give you correct answer even with small number of N
clear;
clc;
close all
f=@(x)2*exp(x).*(2*sin(pi*x)-2*pi*cos(pi*x)+pi^2*sin(pi*x));
sol=@(x)2*exp(x).*sin(pi*x);%exact solution
a=0;
b=1;
N=10;
h=(b-a)/(N+1);
x=(a:h:b)';
x(1)=[];
x(end)=[];
%define diffusion matrix
Ad=2*diag(ones(N,1));
Ad=Ad-diag(ones(N-1,1),-1);
Ad=Ad-diag(ones(N-1,1),1);
Ad=Ad/h^2;
sigma=3;
%define transport matrix
% At=zeros(N,N);
% At=At+diag(ones(N-1,1),1);
% At=At-diag(ones(N-1,1),-1);
% At=At*sigma/(2*h);
At = eye(N,N);
At=At*sigma;
A=Ad+At;
F=f(x);
U=A\F;
U=[0 ; U ; 0];
x=[a ; x ; b];
figure
plot(x,U,'--');%plot approximation
hold on;
plot(x,sol(x));%plot exact solution
legend('app','exact')
Best regards,
Trung
  3 Comments
infinity
infinity on 20 Jun 2019
You are welcome!
I thought that maybe you were confusing between 3u(x) and 3u'(x).
Andrea Cassotti
Andrea Cassotti on 21 Jun 2019
Yes, that is exactly what happened.
Have a nice day!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!