Finite Difference Coding Mistake

13 ビュー (過去 30 日間)
Andrea Cassotti
Andrea Cassotti 2019 年 6 月 20 日
コメント済み: Andrea Cassotti 2019 年 6 月 21 日
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 日
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 件のコメント
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 件)

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by