ode15s solver produces wrong solution

1 view (last 30 days)
Romio
Romio on 29 Jul 2022
Commented: Torsten on 1 Aug 2022
I am trying to solve the heat equation
such that
using the method of line. I discritzed in space and used ode15s to integrate the semi-discrete system of ODEs using ode15s in time, but still when I compare to the exact solution it is wrong. Is there anything wrong in my code?
t0 = 0; % initial time
Tf = 1; % final time
xl = 0; % left point
xr = 1; % right point
m = 52; % no. of ODEs in semidiscrete system
h = 1/m; % space stepsize
k = 0.0001; % time stepsize
Utrue = @(t,x) exp(-t/10) * sin(pi*x); % true solution
xsol = linspace(xl,xr,xr/h);
tspan = linspace(t0,Tf,Tf/k+2);
% initial conditions
U0=zeros(m,1);
for i=1:m
x=i*h;
U0(i) = Utrue(0,x);
end
% solve
opts = odeset('RelTol',1e-13, 'AbsTol',1e-13*ones(m,1),'InitialStep',k, 'MaxStep',k);
[t,u]=ode15s(@ODE,tspan,U0,opts);
% plot solution at x = 0.5 >>> Wrong Solution!!
figure(1)
plot(t,u(:,m/2))
hold on
utrue = Utrue(t,1/2);
plot(t,utrue)
% supporting function
function Ut = ODE(t,U)
m = 52;
h = 1/m; % space stepsize
% BCs
g = zeros(m,1);
Utrue = @(t,x) exp(-t/10) * sin(pi*x); % true solution u
G = @(t,x) -0.1*exp(-t/10)*sin(pi*x) + exp(-t/10)*pi^2 * sin(pi*x); % G = ut-uxx
g(1) = Utrue(t,0)/h^2 + G(t,0);
g(m) = Utrue(t,1)/h^2 + G(t,1);
%semidiscrete system
Ut = zeros(m,1);
K = 1;
Ut(1) = (K/h^2)*(-2*U(1) + 1*U(2));
for i = 2:m-1
Ut(i) = (K/h^2)*(U(i-1) - 2*U(i) + U(i+1));
end
Ut(m) = (K/h^2)*(U(m-1) - 2*U(m));
Ut = Ut + g;
end

Accepted Answer

Torsten
Torsten on 29 Jul 2022
Edited: Torsten on 30 Jul 2022
t0 = 0; % initial time
tf = 1; % final time
nt = 20; % number of output times
xl = 0; % left point
xr = 1; % right point
nx = 51; % no. of ODEs in semidiscrete system
dx = 1/(nx-1); % space stepsize
Utrue = @(t,x) exp(-t/10).*sin(pi*x); % true solution
x = linspace(xl,xr,nx);
tspan = linspace(t0,tf,nt);
% initial conditions
U0 = Utrue(0,x);
% solve
%opts = odeset('RelTol',1e-13, 'AbsTol',1e-14*ones(m,1));
[t,u]=ode15s(@(t,y)ODE(t,y,nx,dx,x),tspan,U0);%,opts);
% plot solution at x = 0.5
figure(1)
plot(t,u(:,(nx+1)/2))
hold on
utrue = Utrue(t,1/2);
plot(t,utrue)
% supporting function
function Ut = ODE(t,U,nx,dx,x)
G = @(t,x) -0.1*exp(-t/10).*sin(pi*x) + exp(-t/10).*pi^2.*sin(pi*x); % G = ut-uxx
%semidiscrete system
K = 1;
Ut = zeros(nx,1);
%Ut(1) = K*(-2*U(1) + 1*U(2))/h^2 + G(t,0);
Ut(2:nx-1) = K*(U(1:nx-2) - 2*U(2:nx-1) + U(3:nx)) / dx^2 + G(t,(x(2:nx-1)).');
%Ut(nx) = K*(U(nx-1) - 2*U(nx))/dx^2 + G(t,1);
end
  2 Comments
Torsten
Torsten on 1 Aug 2022
I know for this problem it is u(t,0) = u(t,1) = 0, but what if it is not?
Then you have a different Utrue (and a different G).
If the boundary conditions at x=0 and x=1 would change with time or whatever, your function "ODE" would be
function Ut = ODE(t,U,nx,dx,x)
Utrue = @(t,x) exp(-t/10).*sin(pi*x); % true solution
G = @(t,x) -0.1*exp(-t/10).*sin(pi*x) + exp(-t/10).*pi^2.*sin(pi*x); % G = ut-uxx
%semidiscrete system
U(1) = Utrue(x(1),t);
U(nx) = Utrue(x(nx),t);
K = 1;
Ut = zeros(nx,1);
%Ut(1) = K*(-2*U(1) + 1*U(2))/h^2 + G(t,0);
Ut(2:nx-1) = K*(U(1:nx-2) - 2*U(2:nx-1) + U(3:nx)) / dx^2 + G(t,(x(2:nx-1)).');
%Ut(nx) = K*(U(nx-1) - 2*U(nx))/dx^2 + G(t,1);
end
Of course, this function also works in your present case, but is not necessary because Ut(1) = Ut(nx) = 0 is set and U(1) = U(nx) = 0 is already coming from the initial conditions for U.

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by