Followup to my previous question. The first order scheme proved unstable for my pde:
$$f_t + A y f_x - B x f_y =0$$
So I'm looking to implement a higher order scheme (using these tables). I was thinking 6th order in time and 4th in space (I'm guessing what the table gives as "Accuracy" is order?). But I don't know how to deal with the boundary conditions? I'm solving the pde on a grid with zero on the boundaries (to model the distribution vanishing in infinity), so what do I do with terms like $f(x-3h,y,t)$ at points where $x=h,2h$? Some sort of extrapolation, or another method altogether, as suggested by David in my previous question?
Here is the current MATLAB code (big thanks to David once again!)
%// Liouville equation clear; %// Equation Parameters: Xmin = -10.0; %// Minimum X Ymin = -10.0; %// Minimum Y Xmax = 10.0; %// Maximum X Ymax = 10.0; %// Maximum Y Tmin=0; %// Minimum time Tmax = 2.5; %// Maximum time A = 1.0; %// A parameter B = 2.0; %// B parameter %// Simulation parameters: Nt = 5000; %// Number of time steps Nx = 150; %// Number of X space steps Ny = 150; %// Number of X space steps %// Create time and space vectors with exactly Nt, Nx, Ny elements t=linspace(Tmin,Tmax,Nt); x=linspace(Xmin,Xmax,Nx); y=linspace(Ymin,Ymax,Ny); %// Make a grid of x and y coordinates [X,Y] = ndgrid(x,y); %// Calculate dt, dx and dy dt=t(2)-t(1); dx=x(2)-x(1); dy=y(2)-y(1); dtdx = dt/dx; %// For simplicity dtdy = dt/dy; u=zeros(Nx,Ny,Nt); %// Initialise u matrix %// Initial condition u(:,:,1)=reshape(mvnpdf([X(:),Y(:)],[1.0 ,0.0],[0.25,0.25]),Nx,Ny); %// Time stepping algorithm i=2:Nx-1; %// Interior spatial points j=2:Ny-1; for k=2:Nt u(i,j,k) = u(i,j,k-1) - 0.5*A*dtdx*Y(i,j).*(u(i+1,j,k-1)-u(i-1,j,k-1)) ... + 0.5*B*dtdy*X(i,j).*(u(i,j+1,k-1)-u(i,j-1,k-1)); end %// Animate the solution filename = 'testnew51.gif'; for m=1:25:Nt contourf(x,y,u(:,:,m)) title(num2str(t(m))) drawnow() frame = getframe(1) im = frame2im(frame); [imind,cm] = rgb2ind(im,256); if m == 1; imwrite(imind,cm,filename,'gif', 'Loopcount',inf); else imwrite(imind,cm,filename,'gif','WriteMode','append'); end pause(0.05) end