CoCalc -- pde_plots.sagews

Project: pde

Views: ¹⁷⁸

Various interacts for plotting solutions to certain partial differential equations

Instruction: to start an interact (or to start it again), press the button. To change an input, make the change and press enter. To change a function (say f) to a constant function (say 1), enter 1+x-x into the f box and press enter. Some time may elapse, so be patient.

The 1-d global heat equation:

$u_t(x,t)-D\,u_{xx}(x,t)=f(x,t)$ , $(x,t) \in \mathbb{R}\times (0,\infty)$ , with initial condition $u(x,0)=g(x)$ for $x\in\mathbb{R}$ .

Duhamel says the solution can be obtained by integration:

u(x,t)=\displaystyle\int_{\mathbb{R}}\Gamma_D(x-y,t)\,g(y)\,dy + \int_0^t\,\int_{\mathbb{R}}\Gamma_D(x-y,t-s)\,f(y,s)\,dy\,ds,

where

\Gamma_D(x,t)=\displaystyle\frac{1}{\sqrt{4\,\pi\,D\,t}}\exp(-\frac{x^2}{4\,D\,t})

is the fundamental solution to

u_t(x,t)-D\,u_{xx}(x,t)=0

(x,t)\in \mathbb{R}\times (0,\infty)

What does the graph of $u$ look like? This depends on the functions $f(x,t)$ , $g(x)$ (which are assumed 0 when $x \not\in [0,1]$ ), and the parameter $D$ . We can use numerical integration methods (Gauss rules, Romberg table) to compute the values of $u(x,t)$ . The interact below shows the graph of $r(x) = u(x,tmax)$ over the interval $\text{[xmin,xmax]}$ or of $q(t)=u(0,t)$ over the interval $\text{[tmin,tmax]}$ . Play with this interact by changing any of the values of f, g, xmin, xmax, tmin, tmax, and D.

y,s,t=var('y,s,t')

@interact(layout={'top': [['D','f','g','c'],['xmin','xmax','tmin','tmax']]})
def getD(D=1/100,f=x+t,g=(x-1/2)^2+1,xmax=3/2,xmin=-1/2,tmax=1/10,tmin=1/1000,c=selector(['t constant', 'x constant','3d'])):
    p=pi.n()
    def romberg_step(R0,a,h,er,f):
        n=len(R0)
        R1=[.5*R0[0]+.5*h*sum([f(a+(2*i+1)*.5*h) for i in range(2^(n-1))])]
        for m in range(1,n+1):
            R1+=[R1[m-1]+1/(4^m-1)*(R1[m-1]-R0[m-1])]
        er=(R1[-1]-R1[-2])
        return R1,er

    def romberg(f,a,b,ER):
        h=b-a
        R0=[.5*h*(f(a)+f(b))]
        k=0
        er=ER
        while k <10 and abs(er)>= ER:
            R0,er=romberg_step(R0,a,h,er,f)
            h=h/2
            k+=1
        return R0,er

    def g1(x):
        return g(x=x)

    def f1(x,t):
        return f(x=x,t=t)
    def GD(x,t):
        return 1/sqrt(4*D*p*t)*exp(-x^2/(4*D*t))
    def u1(x,t):
        return numerical_integral(GD(x-y,t)*g1(y),0,1)[0]

    def u2(x,t):
        def w(s):
            def h(y):
                return GD(x-y,t-s)*f1(y,s)
            return numerical_integral(h(y),0,1)[0]
        bil=romberg(w,tmin,t-.000001,.00001);
        return bil[0][-1]
    def u(x,t):
        return u1(x,t)+u2(x,t)
    def r(x):
        return u(x,tmax)
    x0=(xmax-xmin)/2

    def q(t):
        return u(0,t)
    G=Graphics()
    if c == 't constant':
        mess='time constant at t = '+str(tmax)
        dx=.01*(xmax-xmin)
        pts=([(xmin+i*dx,r(xmin+i*dx)) for i in range(1,99)])
        p=spline(pts)
        G+=plot(p,(xmin,xmax))
        t0=tmax
        tab=[['x'],['u(x,'+str(tmax.n(digits=4))+')']]
        pts2=[]
        x0=xmin
        dx=(xmax-xmin)/4
        for j in range(5):
            pts2+=[pts[j*24]]
            tab[0]+=[pts2[-1][0].n(digits=6)]
            tab[1]+=[pts2[-1][1].n(digits=6)]
            x0+=dx
        G+=point(pts2)
        mess=table(tab)

    elif c=='x constant':
        mess='location constant at x = '+str(x0)
        dt=.01*(tmax-tmin)
        pts=[(tmin+i*dt,q(tmin+i*dt)) for i in range(1,99)]
        p=spline(pts)
        G+=plot(p,(tmin,tmax-.00001))+points(pts)

    else:
        mess='3d plot of u(x,t)'
        G+=plot3d(u,(x,xmin,xmax),(t,tmin+.00001,tmax-.00001),color='grey',opacity=.6,alpha=.5)
    show(G)
    print mess