Search code examples
matlabodenumerical-integration

MATLAB using ode45 to animate wave function throws tolerance error

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.

Solution

  • 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

    enter image description here

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