ode15s solver produces wrong solution

Question

Romio 2022-7-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1770705-ode15s-solver-produces-wrong-solution

评论： Torsten 2022-8-1

采纳的回答： Torsten

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2022-7-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1770705-ode15s-solver-produces-wrong-solution#answer_1017875

编辑：Torsten 2022-7-30

在 MATLAB Online 中打开

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 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 2022-8-1

在 MATLAB Online 中打开

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.

请先登录，再进行评论。