Clear Filters
Clear Filters

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)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!