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

 採用された回答

infinity
infinity 2019 年 6 月 20 日

0 投票

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 件のコメント

Andrea Cassotti
Andrea Cassotti 2019 年 6 月 20 日
Thank you so much! I made such a stupid mistake...
infinity
infinity 2019 年 6 月 20 日
You are welcome!
I thought that maybe you were confusing between 3u(x) and 3u'(x).
Andrea Cassotti
Andrea Cassotti 2019 年 6 月 21 日
Yes, that is exactly what happened.
Have a nice day!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

ヘルプ センター および File ExchangeProgramming についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by