solving the Korteweg-de Vries equation in Matlab

7 visualizaciones (últimos 30 días)
Yoav Eshel
Yoav Eshel el 12 de Dic. de 2020
Comentada: Rahul el 17 de En. de 2025
Hey,
I'm trying to solve the Korteweg-de Vries - in matlab using ode45 and fft, but my code is very inefficient and I figure out why. The odefunc is
function dydt = kdv(y,k)
% du/dt = -d3u/dx3 - 6u du/dx
% numerically approximates the derivative in the fourier domain
yhat = fft(y);
dyhat = 1i*k.*yhat;
d3yhat = -1i*(k.^3).*yhat;
dy = ifft(dyhat);
d3y = ifft(d3yhat);
dydt = - d3y - 6*(y.*dy);
end
and then I call it using
[t,U] = ode45(@(t,y)kdv(y,k),tspan, y0);
but even for very small arrays it takes forever to run. I tried timing each step becuase I suspected that evaluating the fourier transform every iteration is the reason for the inefficiency, but it turns out that this part of the code runs quite fast and the bulk of the computtion time is in ode45 computations between steps.
Any suggestions are appriciated!
  3 comentarios
Yoav Eshel
Yoav Eshel el 20 de Dic. de 2020
I'm trying to solve it without the partial differential equation toolbox. Doesn't the use of fft turns it into an ode?
Rahul
Rahul el 17 de En. de 2025
You are right. The FFT turns it into a ordinary differentuial equation..
the code below does not work , i suspect the problem is in the conv() part
L = 100;
N = 1000;
dx = L/N;
x = -L/2:dx:L/2-dx;
kappa = (2*pi/L)*(-N/2 : N/2-1);
kappa = fftshift(kappa);
u0 = sech(x);
t = 0:0.1:10;
[t,uhat] = ode45(@(t,uhat)kortewegdevries(t,uhat,kappa),t,fft(u0));
for k = 1:length(t)
u(k,:)= ifft(uhat(k,:));
end
figure, waterfall( (u(1:10:end,:)));
figure, imagesc( flipud(u));
function duhatdt = kortewegdevries(t,uhat,kappa)
duhatdt = -1i*(kappa.^3)'.*uhat + conv(uhat,1i*kappa.*uhat); % -(kappa.^2)'.*uhat;%
end

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by