>> N = 2^10; >> dx = 2*pi/N; >> x = -0:dx:2*pi-dx; >> eps = 1; >> u = sqrt(sin(x).^2 + eps); >> plot(x,u), figure(1) >> hold on >> eps = .1; >> u = sqrt(sin(x).^2 + eps); >> plot(x,u), figure(1) >> eps = .01; >> u = sqrt(sin(x).^2 + eps); >> plot(x,u), figure(1) % as epsilon -> 0, the profile develops corners at x=0,pi. Now look % at how this affects the power spectrum >> eps = 1; >> u = sqrt(sin(x).^2 + eps); >> [k,amp] = find_spec(u); >> hold off >> plot(k,amp); figure(1) >> eps = .1; >> u = sqrt(sin(x).^2 + eps); >> [k,amp] = find_spec(u); >> plot(k,amp); figure(1) >> eps = .01; >> u = sqrt(sin(x).^2 + eps); >> [k,amp] = find_spec(u); >> plot(k,amp); figure(1) >> help heat4 [u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N) this solves the periodic heat equation u_t = u_xx with initial data u_0 = cos(x) with periodic boundary conditions using spectral derivatives in space and explicit time-stepping. t_0 is the initial time, t_f is the final time, N is the number of mesh-points, and M is the number of time steps. err is the error. it also returns the power spectrum of the solution, where kk is the collection of wave numbers and amp = log(f_k^2) is the magnitude of the kth Fourier mode. N SHOULD BE A POWER OF 2! (It will work no matter what, but will go much faster if N is a power of 2.) % now we solve the heat equation with explicit time-stepping, having % chosen a time step that won't to lead to trouble. >> [u,err,x,t,kk,amp] = heat4(0,1,100,8); r = 0.0162 r = 0.9100 % plotting the power spectrum, we see that noise stay at the level of % round-off >> plot(kk,amp(:,1)); figure(1) >> hold off; >> plot(kk,amp(:,1)); figure(1) >> hold on >> plot(kk,amp(:,21)); figure(1) >> plot(kk,amp(:,61)); figure(1) >> plot(kk,amp(:,101)); figure(1) % now we solve the heat equation with explicit time-stepping, having % chosen a time step that is sure to lead to trouble. (did this by % taking dx much smaller) >> [u,err,x,t,kk,amp] = heat4(0,1,100,128); r = 4.1501 r = -38.6900 >> clf; hold on; >> plot(kk,amp(:,1)); figure(1) >> plot(kk,amp(:,2)); figure(1) >> plot(kk,amp(:,3)); figure(1) >> plot(kk,amp(:,4)); figure(1) % from the power spectrum, we can see the noise is growing. now we look % for this in the error. >> plot(x,err(:,1)) >> figure(1) >> plot(x,err(:,2)),figure(1) >> plot(x,err(:,3)),figure(1) >> plot(x,err(:,4)),figure(1) >> plot(x,err(:,5)),figure(1) >> plot(x,err(:,6)),figure(1) >> plot(x,err(:,10)),figure(1) >> plot(x,err(:,6)),figure(1) >> plot(x,err(:,7)),figure(1) % up until now, the error was around 1e-5. Now we see that the error % has grown hugely >> plot(x,err(:,8)),figure(1) % looking at the power spectrum, we see that the high modes have grown % to be on the order of 1.e-5, which is why they're now visible in the % error. (They were always present, not just visible.) >> plot(kk,amp(:,8)); figure(1) % now solve with the same dt and dx, but with the implicit method: >> [u,err,x,t,kk,amp] = heat5(0,1,100,128); r = 4.1501 % at this time the power spectrum is still fine. >> plot(kk,amp(:,8)); figure(1) % as is the error: >> plot(x,err(:,8)),figure(1) % now look at the power-spectrum for burger's equation. >> load data.mat >> who Your variables are: N amp_512 k kk_512 u_256 x_256 amp dx kk t u_512 x_512 amp_128 eps kk_128 u x amp_256 err kk_256 u_128 x_128 % plot the solution and watch it become singular in finite time: >> plot(x_512,u_512(:,1)); figure(1); hold on; >> plot(x_512,u_512(:,21)); figure(1); >> plot(x_512,u_512(:,41)); figure(1); >> plot(x_512,u_512(:,61)); figure(1); >> plot(x_512,u_512(:,81)); figure(1); >> plot(x_512,u_512(:,101)); figure(1); >> hold off % plot the power spectrum and watch it expand to fill all modes. >> plot(kk_512,amp_512(:,1)); figure(1); >> hold on >> plot(kk_512,amp_512(:,11)); figure(1); >> plot(kk_512,amp_512(:,21)); figure(1); >> plot(kk_512,amp_512(:,41)); figure(1); >> plot(kk_512,amp_512(:,61)); figure(1); >> plot(kk_512,amp_512(:,81)); figure(1); >> diary off