Implementing numerical method for PDE
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
schlang
am 13 Okt. 2022
Bearbeitet: Davide Masiello
am 13 Okt. 2022
Hello
I am trying to solve the following PDE
with intital and boundary conditions such that
. I used the second order centered finite difference discrtization for
and then want solve the ode system using ode15s
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1154698/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1154703/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1154708/image.png)
Here is my attempt. When I plot the solution obtained from ode15s and compare it to the exact solution they are different. I am not if I made a mistake somewhere. Help is really appreciated
clc,clear,close all
% parameters
t0 = 0;
T = 1.0;
tspan = [t0 T];
xl = 0;
xr = 1;
m = 20;
x = linspace(xl,xr,m + 1);
dx = 1/m;
Uexact = @(t,x) exp(1i*(x-t));
% initial conditions
U0 = Uexact(0,x)';
U0 = U0(2:end-1);
% solve
fn = @(t,U) ODE(t,U,m,dx);
opts = odeset('RelTol',1e-13, 'AbsTol',1e-15);
[t,U] = ode15s(fn, tspan, U0, opts);
%compare with exact solution
plot(x(2:end-1),U(end,:))
hold on
plot(x(2:end-1),Uexact(T,x(2:end-1)))
function dUdt = ODE(t,U,m,dx)
A = eye(m-1);
A = A * (-2);
A = A + diag(ones(m-2,1),1);
A = A + diag(ones(m-2,1),-1);
A = (1/dx^2) * A;
g = zeros(m-1,1);
g(1) = g(1) + (1/dx^2) * exp(1i*(-1*t));
g(end) = g(end) + (1/dx^2) * exp(1i*(1-t));
dUdt = (1i) * (A*U) + g;
end
Thanks
2 Kommentare
Akzeptierte Antwort
Davide Masiello
am 13 Okt. 2022
Bearbeitet: Davide Masiello
am 13 Okt. 2022
clear,clc
tspan = [0 1];
N = 100;
x = linspace(0,1,N);
dx = 1/(N-1);
Uexact = @(t,x) exp(1i*(x-t));
U0 = Uexact(0,x);
M = eye(N);
M(1,1) = 0;
M(N,N) = 0;
opts = odeset('Mass',M,'MassSingular','yes');
[t,U] = ode15s(@(t,U)yourPDE(t,U,N,dx), tspan, U0, opts);
plot(x,real(U(end,:)),'k',x(1:4:end),real(Uexact(1,x(1:4:end))),'r.')
xlabel('x')
ylabel('U')
title('At final time')
legend('Numerical','Exact','Location','Best')
plot(x,real(U(1:3:end,:)),'k',x(1:3:end),real(Uexact(t(1:3:end),x(1:3:end))),'r.')
xlabel('x')
ylabel('U')
title('At several times')
function dUdt = yourPDE(t,U,N,dx)
dUdt(1,1) = U(1)-exp(-1i*t);
dUdt(2:N-1,1) = 1i*(U(1:end-2)-2*U(2:end-1)+U(3:end))/dx^2;
dUdt(N,1) = U(end)-exp(1i*(1-t));
end
0 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Boundary Conditions finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!