

In [3]:

    
#Requires Gadfly
using ApproxFun, Interact

$$u'' + 2u = \cos(\omega t)\qquad u(0)=u'(0)=0$$

Slide to $\omega=\sqrt 2$ to see a resonance, where the solution grows



In [4]:

    
t=Fun(identity,[0.,1000.]) 
S=space(t)
D=Derivative(S)           
B=[ldirichlet(),lneumann()]  

L=D^2+2I  # our differential operator

@manipulate for ω=0.:.1sqrt(2):2sqrt(2)
    u=[B;L]\[0.;0.;cos(ω*t)]
    ApproxFun.plot(u)
end









    











    Out[4]:

Bessel equation $$x^2 u'' + x u + (x^2 - \nu^2)u$$



In [6]:

    
x=Fun(identity,[1.,2000.])  
d=domain(x)                    
D=Derivative(d)               
B=dirichlet(d)  
@manipulate for ν=0:.1:2000
	L=x^2*D^2 + x*D + (x^2 - ν^2)   # our differential operator
   	u=[B;L]\[besselj(ν,first(d)),besselj(ν,last(d))]
   	ApproxFun.plot(u;axis=[-.2,.2])
end









    











    Out[6]:

Exponentially ill-conditioned Lee & Greengard BVP $\epsilon u'' - x u' + u = 0$



In [7]:

    
x=Fun(identity)
d=domain(x)  #[-1.,1.]
D=Derivative(d)
B=dirichlet(d)

@manipulate for ϵ=0.0005:0.0005:0.4
    L=ϵ*D^2 - x*D + I

     u=[B;L]\[1.,2.]
    ApproxFun.plot(u;axis=[-2.,3.])
end









    











    Out[7]:

Piecewise $$u'' + \omega u \begin{cases} 1 & |x| < c\cr 0 & {\rm otherwise}\end{cases} = 0,\qquad u(-1)=1,u(1)=0$$



In [9]:

    
@manipulate for ω=1:50000,c=0.1:0.01:1.
    f=Fun(x->abs(x)<c?1:0,[-1.,-c,c,1.])
    S=space(f)
    B=dirichlet(S)
    D=Derivative(S)    
    ApproxFun.plot([B;D^2+ω*f]\[1.,0.];axis=[-0.5,1.1])
end

Forced Helmholtz equation $$u_{xx} + u_{yy} + k u = {\rm e}^{-(x-x_0)^2-(y-y_0)^2},\qquad \partial u = 0$$



In [10]:

    
d=Interval()^2
B=dirichlet(d)
Δ=lap(d)

@manipulate for ny=10:200,x0=-1.:0.1:1.,y0=-1.:0.1:1.,k=-1000.:1000.
    f=Fun((x,y)->exp(-(x-x0)^2-(y-y0)^2),d)
    u=pdesolve([B;Δ+k*I],[zeros(∂(d));f],ny)
    ApproxFun.contour(u)
end



In [12]:

    
d=Interval()^2
B=[dirichlet(d[1])⊗I,I⊗ldirichlet(d[2]),I⊗rneumann(d[2])]
Δ=lap(d)

@manipulate for ny=10:200,x0=-1.:0.1:1.,y0=-1.:0.1:1.,k=-1000.:1000.
    f=Fun((x,y)->exp(-(x-x0)^2-(y-y0)^2),d)
    u=pdesolve([B;Δ+k*I],[zeros(∂(d));f],ny)
    ApproxFun.contour(u)
end

Dirichlet Helmholtz equation $$\Delta u + k u = 0, u(\pm 1,y)=u(x,\pm 1) = 1$$



In [13]:

    
d=Interval()^2
B=dirichlet(d)
Δ=lap(d)

@manipulate for k=-500.0:.001:2000.0,ny=10:200
    ApproxFun.contour(pdesolve([B;Δ+k*I],ones(∂(d)),ny))
end

Convection diffusion $u_t = \epsilon u_{xx} + (B+C x) u_x, u(x,0)=e^{-20x^2}, u(-1,t)=u(1,t) = 0$



In [14]:

    
dx=Interval();dt=Interval(0,1.)
d=dx*dt
Dx=Derivative(d,1);Dt=Derivative(d,2)
x=Fun(identity,dx)
@manipulate for ε=0.001:0.001:2.,B=-5.:0.1:5.,C=-5.:0.1:5.
    V=B+C*x    
    # Parentheses are a hack to get rank 2 PDE
    u=[timedirichlet(d);Dt-ε*Dx^2-V*Dx]\Fun(x->exp(-20x^2),dx)
    ApproxFun.contour(u;axis=[-1,1,0,1])
end

Schrodinger equation $$i \epsilon u_t = -{\epsilon^2 \over 2} u_{xx} + x^2 u$$



In [15]:

    
dx=Interval(0.,1.);dt=Interval(0.0,0.54)
d=dx*dt

V=Fun(x->x^2,dx)

Dt=diff(d,2);Dx=diff(d,1)

@manipulate for ϵ=0.005:0.002:0.3,ny=50:200
    u0=Fun(x->exp(-25*(x-.5)^2)*exp(-1.im/(5*ϵ)*log(2cosh(5*(x-.5)))),dx)
    L=1im*ϵ*Dt+.5*ϵ^2*Dx^2-V⊗1

    u=pdesolve([timedirichlet(d);L],u0,ny)
    ApproxFun.contour(real(u))
end

Nonlinear BVP $$\epsilon u'' + 6(1-x^2)u' +u^2=1$$ $$u(-1)=u(1)=0$$



In [17]:

    
x=Fun()
u0=0.x
@manipulate for ε=0.001:0.001:1.,c=0.:0.1:1.
    N=u->[u(-1.)-c;u(1.);ε*u''+6*(1-x^2)*u'+u^2-1.]

    u=newton(N,u0)
    ApproxFun.plot(u;axis=[-1.,2.])
end