ode15s solver produces wrong solution

2 次查看（过去 30 天）
Romio
Romio 2022-7-29
评论： Torsten 2022-8-1
Ran in:
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

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-7-29
编辑：Torsten 2022-7-30
Ran in:
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 个评论
Romio
Romio 2022-8-1
@Torsten Thank you very much. I am not sure though how the boundary conditions at x = 0 and x =1 are treated? I know for this problem it is u(t,0) = u(t,1) = 0, but what if it is not?
Torsten
Torsten 2022-8-1
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by