Finite Difference Coding Mistake

1 vue (au cours des 30 derniers jours)
Andrea Cassotti
Andrea Cassotti le 20 Juin 2019
Commenté : Andrea Cassotti le 21 Juin 2019
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

Réponse acceptée

infinity
infinity le 20 Juin 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 commentaires
infinity
infinity le 20 Juin 2019
You are welcome!
I thought that maybe you were confusing between 3u(x) and 3u'(x).
Andrea Cassotti
Andrea Cassotti le 21 Juin 2019
Yes, that is exactly what happened.
Have a nice day!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by