MATLAB using ode45 to animate wave function throws tolerance error

Question

I am working on animating a wave function in MATLAB. So far, I have not reached to part of the animation yet since my ode45 throws an error. I am using de Korteweg-de Vries equations for this (which can be found here: https://en.wikipedia.org/wiki/Korteweg%E2%80%93de_Vries_equation).

Now my code is as follows:

function[t, N, u] = wave(u0, L, T, S, N)

%First we create a vector t which stores the time entries, a vector x which
%stores the location entries and an h which is the location step size for
%the Korteweg-De Vries function.
time = linspace(0,T,S);
x = linspace(-L,L,N);
h = x(2)-x(1);

U0=u0(x);

options=odeset('RelTol',1e-13,'AbsTol',1e-13);
[t,u] = ode45(@kdv,time,U0,options);

    function dudt = kdv(t,u)

    d2udx2 = ([u(N)-2*u(1)+u(2); diff(u,2); u(N-1)-2*u(N)+u(1)] ./ h^2);
    total = d2udx2 + 3.*u.^2;

    diff_total = diff(total);
    diff_total = [diff_total(end); diff_total(1:end)];

    dudt = -[diff_total(2)-diff_total(N-1); diff(diff_total,2); diff_total(N-1)+diff_total(2)] ./ (2*h);

    end
end

Now when I set f=@(x)(2/2).*(1./cosh(sqrt(2).*x./2).^2) and then call the function with [t,N,u]=wave(f, 20, 10, 200, 200); I get the following error:

Warning: Failure at t=6.520003e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.220446e-16) at time t. In ode45 (line 360) In wave (line 23)

Also the returned t and u have a size of 16x1 and 16x200 respectively, but I think they would expand to 200x1 and 200x200 if the warning did not occur.

Any hints on how to solve this would be much appreciated.

What is the reason to increase the iterated differences to order 5 with diff_total? Why is there a second order difference in the final derivative vector? Why do you not make the iterated differences into proper difference quotients by dividing with the appropriate power of h? — Lutz Lehmann
The Korteweg-de vries equations states du/dt = -d/dx * (d^2u/dx^2 + 3 u^2). I therefore think I need the second order difference, to get that part between the brackets, right? — user12545199
Forget my last sentence, the proper divisions are there. d2udx2 and thus total are constructed correctly. Moving additional remarks to answer. — Lutz Lehmann

Lutz Lehmann Lutz Lehmann · Accepted Answer · 2019-12-20T17:04:18

You did too much differences, the middle step of constructing diff_total is redundant.

d2udx2 and thus total are constructed correctly. You could shorten the first construction as
```
d2udx2 = diff([u(N);  u;  u(1) ], 2)./h^2`;
```
The next step constructing diff_total is redundant with the difference construction in dudt. The diff_total construction is lacking a division by h, for some unknown reason you have difference order 2 in the middle of dudt.
The outer x differentiation would be best done via convolution to get a central difference quotient as you have in the outer elements,
```
dudt = -conv([ total(N); total; total(1)], [1; 0; -1], 'valid')/(2*h)
```

All together this gives the derivatives procedure

function dudt = kdv(t,u,h)
    d2udx2 = diff([u(end);  u;  u(1) ], 2)./h^2;
    total = d2udx2 + 3*u.^2;
    dudt = -conv( [ total(end); total; total(1)], [1; 0; -1], 'valid')./(2*h);
end

which (with N=80) produces the solution

which as intended is a soliton wave traveling with speed c=2.

MATLAB using ode45 to animate wave function throws tolerance error

1 Answers