% % finite difference code: % % u + c(x) u = 0 % t x % a = 0.5; % grid refinement, larger "work" means finer grid work = 8; % grids (0 includes zero pt) M = 16*work; h = 2*pi/M; N = 48*work; k = 2*pi/N; x = h:h:2*pi; x0 = [0 x]; t = k:k:2*pi; t0 = [0 t]; % initial value & c(x) vector u = sin(x + a*(1-cos(x))); u0 = [0 u]; c = 1./(1 + a*sin(x)); c0 = [1 c]; fc = (k/h)*c; % solution matrix (begin with IVs) u0_old = 0; u_old = u; u_mn = [0 u]; % iterate on t_n using finite difference formula for n1 = 1:N u0_new = sin(-n1*k); u_new = u_old - fc.*(diff([u0_old u_old])); u_mn = [u_mn; [u0_new u_new]]; u0_old = u0_new; u_old = u_new; end % plot level curves of solution: u_mn = u(x,t) figure(1); clf subplot(2,1,1) plot(x0,u_mn(1,:),'r') hold on plot(x0,u_mn(21:20:161,:),'k') axis([0 2*pi -1.5 1.5]) title('\bf u(x,t) at time intervals of (2\rm\pi\bf/16)') xlabel('\bf x'); ylabel('\bf u(x,t), t=0 is red') subplot(2,1,2) plot(x,fc) axis([0 2*pi 0 2]) title('\bf CFL number (k*c(x)/h)') xlabel('\bf x'); ylabel('\bf CFL # \rm(>1, unstable)') axis([0 2*pi 0 2]); hold on plot([0 9],[1 1],'--')