bipolar + wavelets

Complex/DirichletBipolar.m

@@ -0,0 +1,32 @@
+function [] = DirichletBipolar(R1, R2, L, N)
+N(1:2)=N;
+a=sqrt((R1^2-R2^2)^2-8*L^2*(R1^2+R2^2)+16*L^4)/(4*L);
+x1=-asinh(a/R1);
+x2= asinh(a/R2);
+
+[Dx,x]=chebD(N(1)); x=x1+(x2-x1)/2*(x+1); Dx=2/(x2-x1)*Dx; Dxx=Dx*Dx; 
+[Dyy,y]=periodicD2(N(2));
+[xx,yy]=ndgrid(x,y);
+zz=xx+1i*yy;
+ww=a*coth(zz/2)-a*coth(x2)+L;
+uu=real(ww);
+vv=imag(ww);
+J=abs(-a/2*csch(zz/2).^2).^2;
+F=zeros(N);
+
+bc=[1+0*y; -1+0*y];
+RHS=J.*F-Dxx(:,[1,end])*bc;
+psi=zeros(N);
+psi([1,end],:)=bc;
+psi(2:end-1,:)=sylvester(Dxx(2:end-1,2:end-1), Dyy', RHS(2:end-1,:));
+
+h=mesh(uu(:,[1:end,1]),vv(:,[1:end,1]),psi(:,[1:end,1]));
+colormap(jet(256)); shading interp;
+set(h,'FaceColor','none'); whitebg(1,'k'); view(2);
+c=4*(R1+R2+L); xlim([-c,c]); ylim([-c,c]);
+xrange=2*c;
+yrange=2*c;
+zrange=max(psi(:))-min(psi(:));
+daspect([1 1 2*zrange/hypot(xrange,yrange)]);
+end
+

Complex/bipolarMap.m

@@ -0,0 +1,12 @@
+function [ww] = bipolarMap(R1,R2,L,N)
+N(1:2)=N;
+a=sqrt((R1^2-R2^2)^2-8*L^2*(R1^2+R2^2)+16*L^4)/(4*L);
+x1=asinh(a/R1);
+x2=asinh(a/R2);
+
+x=linspace(-x1,x2,N(1));
+y=linspace(0,2*pi,N(2));
+[xx,yy]=ndgrid(x,y);
+zz=xx+1i*yy;
+ww=a*coth(zz/2)-a*coth(x2)+L;
+end

Complex/ellipticMap.m

@@ -0,0 +1,10 @@
+function [ww] = ellipticMap(a,b,N)
+N(1:2)=N;
+f=sqrt(a^2-b^2);
+xi0=acosh(a/f);
+x=linspace(0,xi0,N(1));
+y=linspace(0,2*pi,N(2));
+[xx,yy]=ndgrid(x,y);
+zz=xx+1i*yy;
+ww=f*cosh(zz);
+end

Image/imhaart.m

@@ -0,0 +1,28 @@
+function [] = imhaart()
+% Black and white image
+X=imread('Lenna.png');
+X=double(X)/255;
+
+% Gaussian noise
+sigma=0.08;
+w=sigma*randn(size(X));
+Y=max(min(X+w,1),0);
+
+% Haar wavelet transform
+N=wmaxlev(min(size(X(:,:,1))),'haar');
+[C,S]=wavedec2(Y,N,'haar');
+
+% Kill lower coefficients
+tol=3*sigma;
+CT=C;
+CT(abs(C)<tol)=0;
+CT(1:3*S(1,1)*S(1,2))=C(1:3*S(1,1)*S(1,2));
+
+% Inverse Haar wavelet transform
+Z=waverec2(CT,S,'haar');
+
+% Show images
+figure(1); 
+imagesc([X,Y,Z]);
+axis equal off;
+end

ODE/nonlinearODE.m

@@ -1,25 +1,20 @@
 function [] = nonlinearODE(n)
 [D,x]=chebD(n);
 D2=D*D;
-
 u1=1;
 u2=2;
 
 bc=D2(2:end-1,[1 end])*[u2;u1];
 u=u1+(u2-u1)*(x+1)/2;
 
-
-
-i=0;
-err=1;
-tol=1e-10;
+i=0; err=1; tol=1e-10;
 while err>tol
     unew=D2(2:end-1,2:end-1)\(exp(u(2:end-1))-bc);
     err=norm(unew-u(2:end-1),inf);
     u(2:end-1)=unew;
     i=i+1;
 end
 fprintf('%d iterations\n',i);
-h=plot(x,u);
+plot(x,u);
 end
 

Physics/brachistochrone.m

@@ -12,10 +12,16 @@
 opts=[];
 opts.FunctionTolerance=1e-13;
 
-[path,t]=ga(@timeFun,n-2,[],[],[],[],min(y1,y2)-L,max(y1,y2)+L,[],[],opts);
+lb=(min(y1,y2)-L)*ones(n-2,1);
+ub=(max(y1,y2)+L)*ones(n-2,1);
+path=y1+(y2-y1)*(x(2:end-1)+1)/2;
+[path,t]=ga(@timeFun,n-2,[],[],[],[],lb,ub,[],[],opts);
+disp(t);
+[path,t]=fmincon(@timeFun,path,[],[],[],[],lb,ub);
+disp(t);
+
 y=[y2;path(:);y1];
 plot(x,y);
-disp(t);
 end
 
 function t=timeFun(y)

Physics/pwePML1.m

@@ -11,7 +11,7 @@
 % Layer
 xl=3/4*x(end);
 layer=abs(x)>xl;
-roi=1:N;%~layer;
+roi=~layer;
 sig=zeros(N,1);
 sig(layer)=2/dz*((abs(x(layer))-xl)/(x(end)-xl)).^3;
 
@@ -23,7 +23,7 @@
 
 figure(1);
 h=plot(x(roi), w(roi,1));
-xlim([-xl,xl]);
+%xlim([-xl,xl]);
 ylim([-1,1]);
 
 nframes=10000;

Physics/pwePML2.m

@@ -6,6 +6,7 @@
 lambda=100;
 k=2*pi/lambda;
 [D, x]=hermD(N);
+
 [xx,yy]=ndgrid(x);
 dz=6/N^2;
 
@@ -36,9 +37,8 @@
     w(:,[1 end],:)=0;
     if(mod(i,5)==1)
         E=w(roi,roi,1);
-        EE=real(E).^2+imag(E).^2;
 
-        V=mat2gray(EE);
+        V=mat2gray(abs(E).^2);
         V=cat(3,V,V,V);
         phase=angle(E)/(2*pi);
         phase=phase-floor(phase);

Physics/sinegordon.m

@@ -0,0 +1,24 @@
+function [] = sinegordon( N )
+[D,x]=chebD(N); t=1+x';
+Dxx=D*D;
+Dtt=D*D;
+Dtt(1,:)=D(end,:);
+
+u=zeros(N);
+ut0=-0*exp(-10*x.^2);
+u(:,1)=ut0;
+u(:,end)=exp(-10*x.^2);
+BC=Dxx*u-u*Dtt';
+
+i=0; err=1; tol=1e-10;
+while err>tol
+    unew=sylvester(Dxx(2:end-1,2:end-1), -Dtt(1:end-1,1:end-1)', [ut0(2:end-1), sin(u(2:end-1,2:end-1))]-BC(2:end-1,1:end-1));
+    err=norm(unew-u(2:end-1,1:end-1),inf);
+    u(2:end-1,1:end-1)=unew;
+    i=i+1;
+
+end
+[xx,tt]=ndgrid(x,t);
+surf(xx,tt,u); colormap(jet(256)); view(2);
+end
+

Signal/noise.m

@@ -1,7 +1,7 @@
 function [pp] = noise(p,t)
-p=p+randn(size(p));
-p_hat=fft(p);
-pp=sqrt(2*pi)/length(p)*ifft(conj(p_hat).*p_hat);
-plot(t, [pp; cos(t)*sqrt(pi/2)]);
+pr=p+randn(size(p));
+p_hat=fft(pr);
+pp=2/length(p)*ifft(conj(p_hat).*p_hat);
+plot(t, [pp; p]);
 end