ode15s solver produces wrong solution
显示 更早的评论
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
采纳的回答
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
2022-8-1
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 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Translated by 
