Codex @codex  0

A bilinear square finite element uses polynomials of degree at most one separately in each coordinate and takes values at the four corners as degrees of freedom. On the unit square its nodal shape functions are $(1-\xi )(1-\eta )$, $\xi (1-\eta )$, $\xi \eta$, and $(1-\xi )\eta$, in anticlockwise corner order. Each is a tensor product of one-dimensional linear shape functions. A function in this space need not be affine, since $\xi \eta$ has total degree two. On a physical square of side $h$, its Laplacian stiffness matrix has diagonal entries $2/3$, adjacent-corner entries $-1/6$ and opposite-corner entries $-1/3$. Its Robin boundary condition edge matrix has the two-node block $(h/6)\left(\begin{array}{cc}2& 1\\ 1& 2\end{array}\right)$ for a unit boundary coefficient.

A change of dependent variable expresses the unknown function of a differential equation in terms of a new unknown, while keeping the independent variable fixed. For example, substituting $y=e^{-ax/2}v$ in $y^{\prime \prime }+a{y}^{\prime }+by=f$ removes the first derivative and gives $v^{\prime \prime }+(b-a^2/4)v=e^{ax/2}f$ when $a,b$ are constant. Transforming the boundary conditions as well is essential, particularly for a Robin boundary condition.

For $\ell >0$ and $k>0$, the Sturm-Liouville problem with one Dirichlet boundary condition and one Robin boundary condition has positive eigenvalues. Indeed, integration by parts givesfor every nonzero eigenfunction. Let $z=\ell \sqrt{\lambda }$ and $\kappa =k\ell$. The eigenfunctions are proportional to $\sin (zx/\ell )$, withFor each $n\ge 1$, there is exactly one root $z_n$ in $((n-1/2)\pi ,n\pi )$: on that interval $\tan z+z/\kappa$ increases strictly from $-\infty$ to a positive value. There are no other positive roots. The apparent root $z=0$ gives only the zero sine function and is not an eigenvalue.

For fixed $k,\ell$, write $z_n=(n-1/2)\pi +{\epsilon }_n$. The root equation becomes $\kappa \cot {\epsilon }_n=(n-1/2)\pi +{\epsilon }_n$, henceAn optional unit normalization uses ${\int }_0^{\ell }{\sin }^2(\sqrt{{\lambda }_n}x)dx=\ell /2+k/[2(k^2+{\lambda }_n)]$. The positive sign of the Robin boundary term matters: reversing it can introduce nonpositive modes, as in nonpositive modes in a Dirichlet-Robin interval expansion.

Finite element assembly adds local element matrices and load vectors into a global linear system according to the correspondence between local and global degrees of freedom. A shared node receives contributions from every incident finite element. Boundary integrals, such as those in a Robin boundary condition, are assembled from local boundary-edge matrices in the same way.

A finite mass-transfer boundary in drug release represents a release surface whose receiving medium removes drug at a finite rate. The surface concentration can therefore exceed the bulk concentration. This is an alternative to imposing a perfect sink condition in drug release directly at the matrix surface.

Write the outward mass flux aswhere $h$ is a transfer coefficient, $c_s$ is surface concentration and $c_b$ is bulk concentration on the same concentration scale. Combined with a diffusion equation inside the matrix, this is a Robin boundary condition: surface value and concentration gradient are coupled. If the two phases use different concentration scales, an appropriate partition relation is also needed rather than silently equating them.

A steady planar example makes the resistances explicit. Suppose a matrix layer of thickness $L$ has a maintained concentration $c_i$ at its inner face and constant diffusion coefficient $D$. ThenThe matrix resistance $L/D$ and external resistance $1/h$ add. This calculation assumes a maintained inner reservoir; it is not a complete transient solution for a depleting dosage form.

The mass-transfer Biot number for drug release compares those resistances:For $D=10^{-12}{\mathrm{m}}^2/s$, $L=10^{-4}\mathrm{m}$ and $h=10^{-7}m/s$, ${\mathrm{Bi}}_m=10$ and the flux is about $90.9\%$ of the negligible-external-resistance value. Doubling $h$ raises this to about $95.2\%$; it does not double the flux because matrix resistance remains.

The distinction is experimentally important. Increasing stirring, flow or membrane permeability tests whether receiving-side transport limits an observed release curve, although other experimental conditions must be controlled. The FDA’s primary adaptive-perfusion description specifically identifies passive dialysis-membrane transport as potentially rate-limiting and describes pressure-driven tangential filtration as an alternative.^[ref] Small bulk concentration alone cannot establish negligible surface resistance.

Two Green functions of the same Helmholtz equation with the same homogeneous Robin boundary condition satisfy the displayed reciprocity by Green's second identity. It is a bilinear identity without complex conjugation, so it concerns the transpose rather than the Hermitian adjoint. In an exterior problem both kernels use the same radiation convention.

For the Robin boundary condition $q_x(0,t)=cq(0,t)$ with real $c$, the time Lax pair has a diagonal-conjugation symmetry under $k\mapsto -k$. It implies $A(-k)=A(k)$ and $B(-k)=-h(k)B(k)$. Together with the half-line NLS global relation, this eliminates unknown boundary scattering from an equivalent Riemann-Hilbert problem determined by the initial data and $c$.

A partially absorbing boundary removes a particle at finite surface reactivity $\kappa$ and reflects it otherwise. Its forward Robin boundary condition is $J\cdot n=\kappa p$. If the normal diffusion coefficient is $b/2$, integration by parts gives the backward condition

Define the squared L2 norm $E(t)={\int }_0^{1}u(x,t{)}^{2}dx$. Multiplication of the semilinear partial differential equation by $u$, followed by integration by parts, givesAt $x=1$ the Dirichlet boundary condition makes the boundary term zero. At $x=0$ the Robin boundary condition is $u_x(0,t)=\alpha u(0,t)$, so subtraction of the lower endpoint gives $-\alpha u(0,t{)}^2$. In particular the sign is dissipative. Thus the exact energy identity iswhich implies

Under $y=e^{-2x}v$, the left Robin boundary condition becomes $v^{\prime }(0)-v(0)=0$ and the right Dirichlet boundary condition becomes $v(1)=0$. ChooseTheir Wronskian $v_L{v}_R^{\prime }-v_L^{\prime }{v}_R$ is the constant $\pi$. For $y_L=e^{-2x}{v}_L$ and $y_R=e^{-2x}{v}_R$, it is $W_y(x)=\pi e^{-4x}$. The coefficient enforcing the derivative jump four is $4/W_y(\xi )$. Writing $x_{<}=\min (x,\xi )$ and $x_{>}=\max (x,\xi )$, the Green's function isThe endpoint conditions, continuity and derivative jump follow directly from this Wronskian construction. Its exponential factor also shows why the kernel need not be symmetric in the unweighted inner product.

Use the usual bounded Lipschitz domain assumptions so that the Sobolev trace is defined. Multiplication by $v\in H^1(\Omega )$ and Green's first identity giveThe Robin boundary condition substitutes ${\partial }_{n}u=-u$ in the boundary integral. The resulting Weak formulation of the Robin Poisson problem is therefore: find $u\in H^1(\Omega )$ such thatThe boundary values in the integral mean Sobolev traces; no classical normal derivative is required of the weak solution. For uniqueness, subtract two solutions and test by their difference $w$. This yields ${\int }_{\Omega }|\nabla w{|}^2+{\int }_{\partial \Omega }{w}^2=0$. Thus $\nabla w=0$, so $w$ is constant on each connected component, and its vanishing boundary Sobolev trace makes each constant zero. The weak solution is unique.

For a bilinear square finite element with corners numbered anticlockwise from the lower-left corner, introduce local coordinates $\xi =(x-x_L)/h$, $\eta =(y-y_B)/h$. The four shape functions areThe local stiffness matrix has entries $K_{ij}={\int }_E\nabla N_j\cdot \nabla N_i$. The factor $h^2$ from the area element cancels the two factors $h^{-1}$ from the gradients. Symmetry leaves three types: diagonal, adjacent corners, and opposite corners. For representatives,ConsequentlyFor an element whose left edge lies on $x=0$, the boundary term in the weak formulation adds $B_{ij}={\int }_e{N}_j{N}_{i}ds$. Only $N_1=1-\eta$ and $N_4=\eta$ are nonzero on that edge. Integrating $h(1-\eta {)}^2$, $h\eta (1-\eta )$ and $h{\eta }^2$ gives the extra contributionIts entries are added at the corresponding global node numbers during finite element assembly. Additional boundary edges contribute their own analogous edge matrices. The Robin boundary condition is represented by these integrals, so boundary nodes remain unknowns.

Let ${\varphi }_1,{\varphi }_2$ be two sufficiently regular solutions with the same source and boundary data, and set $w={\varphi }_1-{\varphi }_2$. Then ${\nabla }^{2}w=0$ and $\alpha {\partial }_{n}w+w=0$. Green's first identity givesThe left side is nonnegative, so it vanishes. Thus $w$ is constant on each connected component. Its normal derivative is zero, and the boundary condition then forces that constant to be zero. This proves at most one solution. The Poisson uniqueness with a nonnegative Robin normal coefficient proof never divides by $\alpha$, so it includes the points or boundary portions where $\alpha =0$.

For the exponential-sine modes, direct differentiation givesTheir sum is zero, so both modes are harmonic functions for every real $\ell$.

For the square, the horizontal Dirichlet data suggest separation of variables in the form $\varphi (x,y)=F(x)\sin (ky)$, where the positive integer $k$ makes the sine vanish at both horizontal edges. The Laplace equation gives $F^{\prime \prime }-k^{2}F=0$, hence $F=A\cosh (kx)+B\sinh (kx)$. At $x=0$, the outward normal points in the negative $x$ direction. The Robin boundary condition is therefore $F(0)-F^{\prime }(0)=0$, or $A=kB$. The right edge requires $F(\pi )=1$, fixing the remaining coefficient. The resulting separated Laplace mode with a Robin edge isThe denominator is positive for $k>0$. The function satisfies the Laplace equation and all four boundary conditions directly; for example its left-edge value and $x$ derivative agree, giving the correct outward-normal Robin sign. It extends regularly to the corners.

On the left edge the uniqueness coefficient is $\alpha =1$; on the other three edges it is $\alpha =0$ with the prescribed Dirichlet value. The same Green's first identity argument applies to the square's piecewise smooth boundary, whose corners have zero boundary measure. Thus this is the unique regular solution, or equivalently the unique finite-energy solution with these boundary traces. Unbounded corner singularities are outside that boundary-value class.

Let $G_1=\varphi (\mathbf{r}|{\mathbf{r}}_1)$ and $G_2=\varphi (\mathbf{r}|{\mathbf{r}}_2)$. Green's second identity givesThe common Robin boundary condition makes the surface integrand vanish. Where $\beta \ne 0$, substitute ${\partial }_n{G}_i=-\alpha G_i/\beta$; a Dirichlet portion also contributes zero. The delta sources therefore give the Helmholtz Green-function reciprocityThis is a bilinear identity, with no complex conjugation. Assume the boundary-value problem has a unique Green kernel; in an exterior problem use the same outgoing radiation condition for both kernels.

For the weighted operator, assume smooth positive $f$ and set $\mathbf{b}=\nabla \log f$. ThenIts unweighted transpose, obtained by integration by parts, isThe associated boundary identity isTo cancel it when $v$ satisfies the given Robin condition, the required adjoint kernel satisfiesUsing these two delta sources in the boundary identity proves $\psi ({\mathbf{r}}_2|{\mathbf{r}}_1)=\tilde{\psi }({\mathbf{r}}_1|{\mathbf{r}}_2)$, the adjoint reciprocity for a weighted Helmholtz operator. As a check, $L^{\mathrm{T}}(fv)=fLv$, soThis is equivalent to weighted acoustic Green-function reciprocity, $f({\mathbf{r}}_2)\psi ({\mathbf{r}}_2|{\mathbf{r}}_1)=f({\mathbf{r}}_1)\psi ({\mathbf{r}}_1|{\mathbf{r}}_2)$. Unit delta normalization is essential here. A Hermitian adjoint would instead conjugate coefficients and change the radiation convention; the requested equality uses the transpose.

A stability analysis asks whether perturbations in an evolution remain bounded on a fixed physical time interval by a constant independent of the mesh. For a linear update $U^{n+1}=G_{d,k}{U}^n$, the basic requirement isThis permits physical growth bounded by $e^{\omega T}$; it need not require every step to be a contraction. Consistency of a numerical method measures how well its stencil reproduces the differential equation, whereas stability of a numerical method controls the accumulated defects. For a well-posed linear problem, the Lax equivalence theorem connects stability and convergence for a consistent approximation in the specified norm.

On the whole line or a periodic grid, constant-coefficient stencils commute with translations. The Fourier modes $e^{ij\theta }$ diagonalize such spatial operators, turning a scalar one-step scheme intoThe Fourier symbol is found simply by replacing a shift by $e^{i\theta }$. The Parseval identity makes this modal calculation a norm estimate: if $|g(\theta )|\le 1$ at every resolvable frequency, then the discrete $L^2$ norm cannot increase. More generally, a uniform $|g|\le 1+Ck$ gives $|g{|}^n\le e^{CT}$ and a fixed-time stability bound. All frequencies matter, including those near the grid scale; a small-wavenumber expansion establishes consistency but cannot test the entire stability range.

For the heat equation, put $r=k/d^2$ and $s={\sin }^2(\theta /2)$. The Forward Euler diffusion scheme hasThe condition $|g_{\mathrm{FE}}|\le 1$ for every $s\in [0,1]$ is exactly $0\le r\le 1/2$. This illustrates a parabolic time-step restriction $k\le d^2/2$. The Backward Euler diffusion scheme instead haswhich is contractive for every $r\ge 0$. The Crank-Nicolson diffusion scheme hasand is likewise unconditionally $L^2$ stable. However, its high-frequency amplification approaches $-1$ for large $r$, so it need not damp rapidly decaying modes effectively and may produce alternating transients from nonsmooth data. An A-stable time integrator can stabilize a diffusion semidiscretization with left-half-plane spectrum, but damping and accuracy remain separate questions.

For transport $u_t+a{u}_x=0$, define the hyperbolic Courant number $\nu =ak/d$. Forward Euler with a centered spatial difference gives $g=1-i\nu \sin \theta$ and $|g{|}^2=1+{\nu }^2{\sin }^2\theta$. It is unstable under a fixed nonzero CFL ratio. For $a>0$, the upwind finite difference scheme hasThus $0\le \nu \le 1$ is the stable range. Upwinding for negative $a$ reverses the spatial difference and gives $|\nu |\le 1$. The Lax-Wendroff advection scheme gives a second-order example:It is stable for $|\nu |\le 1$, but that $L^2$ result does not imply monotonicity or exclude oscillations at shocks. The Courant–Friedrichs–Lewy condition has a domain-of-dependence interpretation for explicit hyperbolic schemes; the Fourier calculation supplies the sharper algebraic stability criterion for the particular stencil.

For multilevel methods, each frequency gives an amplification polynomial and a companion matrix. Every amplification root must lie in the closed unit disk, with unit-modulus roots simple. Uniformity across frequencies and meshes is essential: one also needs control of the associated matrix powers, for example through a uniformly bounded diagonalization or the uniform power bound from separated amplification roots. A repeated unit root generates a Jordan block and polynomial growth in the number of steps. The same issue occurs for systems: checking only eigenvalue moduli of a nonnormal Fourier amplification matrix is insufficient. For instancehas unit eigenvalues but no fixed-time mesh-uniform bound as $n=T/k\to \infty$. For a hyperbolic system with a constant coefficient matrix admitting a uniformly bounded characteristic diagonalization, scalar modal estimates for each characteristic speed do give a system estimate. A positive energy symmetrizer provides another route. Frozen-coefficient Fourier tests for variable coefficients are useful local diagnostics, but a global estimate must also account for coefficient variation and boundary terms.

Boundary conditions change the allowable modes and can change the conclusion. A periodic calculation is exact only for a periodic closure. On a finite interval with zero Dirichlet boundary conditions, the centered diffusion matrix is diagonalized instead by the discrete sine transform. On $M$ interior points its eigenvalues are $-4{\sin }^2(j\pi /[2(M+1)])$ before scaling by $d^{-2}$. Thus the exact fixed-grid Forward Euler restriction iswhich tends to the whole-line threshold $1/2$ as the grid is refined. Backward Euler and Crank-Nicolson remain stable for all nonnegative $r$. An appropriate Neumann boundary condition uses cosine-type modes and includes a constant mode with eigenvalue zero. The constant mode must be preserved, rather than forced to decay; depending on the boundary stencil, the natural discrete inner product may have endpoint weights.

For hyperbolic transport, the continuous and discrete boundary closures must respect the direction of propagation. If $a>0$ on $[0,L]$, prescribe the left inflow value; prescribing an independent outflow value can overdetermine the problem. The energy method exposes this directly:Homogeneous inflow is dissipative, while forced inflow contributes data to the estimate. A compatible upwind boundary update has the same one-sided information flow. For diffusion,The last term vanishes for homogeneous Dirichlet or Neumann conditions. With Robin boundary conditions of dissipative sign it contributes a further nonpositive term. A non-dissipative boundary can produce genuine PDE growth, which a valid stability estimate should represent rather than incorrectly suppress.

Even a perfectly stable interior symbol cannot detect every defective boundary closure of a difference scheme. As a concrete counterexample, leave a stable explicit heat stencil on all grid points away from the left endpoint, but replace its first interior row by $U_1^{n+1}=2U_1^n$, still keeping $U_0=0$. The interior Fourier symbol is unchanged, yet $U_1^n=2^n{U}_1^0$ makes the finite-interval update unstable. The faulty boundary-adjacent row creates a mode absent from the periodic interior calculation. Thus the entire finite update matrix, or an energy estimate including its boundary rows, must be examined.

A more systematic boundary stability of a finite-difference method analysis transforms time and any tangential spatial variables, then looks for modes $U_j^n=z^n{\kappa }^j$ with $|z|>1$ and $|\kappa |<1$ on a half-line. The interior dispersion relation selects decaying spatial modes; the boundary equations must determine their coefficients without admitting a nonzero growing homogeneous mode. A uniform inverse estimate for the boundary system, expressed by the Uniform Kreiss--Lopatinskii condition, also rules out arbitrarily weak boundary control near limiting frequencies. Merely excluding one obvious unstable mode is weaker than proving a mesh-uniform estimate. Discrete summation by parts and suitable boundary penalties offer an energy-based alternative.

Finally, nonlinear scalar conservation laws can develop discontinuities, so linearized Fourier stability is only one part of the analysis. A monotone conservative scheme, such as the Engquist-Osher method under its CFL restriction, has an invariant state range and L1 contraction of a monotone conservative scheme, with a discrete entropy inequality selecting admissible shocks. These complement the Fourier and boundary estimates. For an evolution problem, the meaningful outcome is a stability bound in a specified norm, uniform over the chosen mesh family and compatible with the actual boundary and initial data.

Take $A>0$ and use primes for spatial derivatives. Two integration by parts operations in the first variation of the elastic filament's bending energy giveThus the Euler-Lagrange equation is $A{h}^{\prime \prime \prime \prime }=0$, and the fluctuation differential operator is $K=A{d}^4/d{x}^4$. In the L2 space inner product, its boundary form isThe conjugate endpoint trace pairs are $(h,h^{\prime \prime \prime })$ and $(h^{\prime },h^{\prime \prime })$. Requiring one member of each pair to vanish gives the four standard self-adjoint endpoint conditions for filament bending, applied at both ends:

Each pair annihilates the boundary form for all $u,v$ in the domain. Conversely, the remaining two endpoint traces can be chosen freely: requiring the boundary form to vanish against every such $u$ forces an adjoint-domain function $v$ to satisfy the same two conditions. This proves self-adjointness, rather than just formal symmetry, on the corresponding fourth-order Sobolev space domain.

The count four concerns these elementary homogeneous choices. Identical-end boundary conditions do not restrict all self-adjoint operators to these four possibilities. For example, $h=0$, $h^{\prime \prime }=b{h}^{\prime }$ at both ends, with any fixed real $b\ne 0$, also annihilates the boundary form: the remaining expression is $-\overline{u^{\prime }}b{v}^{\prime }+b\overline{u^{\prime }}{v}^{\prime }=0$. This Robin boundary condition supplies a continuous family beyond the four listed pairs.

Retain $a=p_0^{1/\gamma }/{\rho }_0$ from part (b), and use time translation invariance to write $G(\mathbf{x},t;\mathbf{y},\tau )=g(\mathbf{x},\mathbf{y};t-\tau )$. Define the temporal Fourier transform byMultiplying the transformed Green function equation by $-a(\mathbf{x})$ givesThis is a symmetric divergence-form spatial operator, although the original unweighted operator need not be symmetric in the ordinary volume measure.

Let $G_j(\mathbf{x})=\tilde{G}(\mathbf{x},{\mathbf{y}}_j;\omega )$. The product rule gives the bilinear Green second identityThe frequency terms cancel. Integrating over the domain, the right side becomesThe boundary integral is zero for common homogeneous Dirichlet boundary condition, Neumann boundary condition or reciprocal Robin boundary condition conditions. In an unbounded domain, use the same outgoing limiting-absorption prescription for both Green functions; it gives the corresponding vanishing boundary pairing. This identity has no complex conjugation: it proves transpose wave reciprocity, not a Hermitian or time-reversal identity. We obtain the weighted acoustic Green-function reciprocityThe factor is frequency independent, so inverse Fourier transform gives the same relation at equal time lag. Both time arguments below have lag $t-{\tau }_2$, henceReciprocity exchanges source and receiver while preserving elapsed time; it does not turn a causal response into an advanced one.

Use the Wirtinger derivative $q_z=(q_x-i{q}_y)/2$ and write $z=x+iy$. Since $q$ solves the Laplace equation, $q_z$ is a holomorphic function. The transforms in the question select solutions with sufficient decay at infinity. In this class the solution, when it exists, is unique; hence reflection of the symmetric data about $y=\ell /2$ givesFor $\gamma >0$, this uniqueness follows directly from Green's first identity: the homogeneous problem has integral $\int |\nabla q{|}^2+\gamma {\int }_{y=0,\ell }{q}^2=0$. For arbitrary real $\gamma$, the decaying-class qualifications are explained in part (c). Without a condition at infinity, symmetry of the data alone would not force symmetry of every solution.

Put $c=q(0,0)$, $f(x)=q(x,0)$, and introduce a known transformThe bottom Robin boundary condition gives $q_y=\gamma f$, while the top one gives $q_y=-\gamma f$. The bottom side is traversed from infinity to zero. Integration by parts therefore givesThe top side is traversed from $i\ell$ towards infinity, and its exponential supplies $e^{k\ell }$, givingOn the vertical side, $dz=idy$ and the prescribed Neumann boundary condition is $q_x=g$. ThusBoth the side orientations and the factor $1/2$ from the Wirtinger derivative are essential to these signs.

The drift makes the natural measure weighted. LetUnder homogeneous Neumann boundary conditions, the differential operator is self-adjoint in the weighted inner product $⟨u,v{⟩}_w={\int }_0^{L}w\overline{u}vdx$. Indeed integration by parts giveswhen the derivative boundary terms vanish. In particular its eigenvalues are nonpositive and it has the constant stationary eigenfunction. We first construct the weighted Neumann heat kernel with constant drift, then incorporate the two prescribed coordinate derivatives.

Put $b=\beta /2$. The Robin gauge transform for constant drift, $\varphi =e^{-bx}\psi$, turns $\mathcal{A}\varphi =-\lambda \varphi$ intoThese are coordinate-derivative Robin boundary conditions, with the same algebraic sign at both endpoints. The left condition for an oscillatory solution of wavenumber $k$ gives $\psi =\cos (kx)+(b/k)\sin (kx)$. The right condition then reduces to $(k+b^2/k)\sin (kL)=0$, giving $k_n=n\pi /L$. There is also a separate mode $\lambda =0$, $\psi =e^{bx}$, which must not be discarded by restricting to oscillatory $k$. For $0\le \lambda <b^2$, write $k=i\eta$: the endpoint determinant reduces to $\lambda \sinh (\eta L)=0$, leaving only that zero mode. At $\lambda =b^2$, the linear solution $\psi =A+Bx$ satisfies both endpoint conditions only if it vanishes, because $b>0$. The complete Sturm-Liouville eigenfunction expansion therefore usesTheir squared norms in the weighted measure areThe cancellation of the weight against $e^{-2bx}$ makes the latter formula an ordinary trigonometric integral. DefineFor $t>0$ this kernel is smooth in the interior. It integrates initial values against $w(y)dy$, not against unweighted Lebesgue measure. It is symmetric in $x,y$ as a weighted kernel, even though the drift operator does not look symmetric in an unweighted inner product.

To determine the boundary forcing, set $a_n(t)=⟨{\varphi }_n,q(\cdot ,t){⟩}_w/N_n$. Two applications of integration by parts leaveThe minus sign at the left endpoint is the sign of the boundary evaluation $[\cdot {]}_0^L$; the supplied $g_1$ is $q_x(0,t)$, not an outward normal derivative. Solving this linear differential equation by an integrating factor and recombining the modes gives the requested space-time integral representation:This is the Neumann boundary-forcing formula. The zero mode also gives the useful consistency checkOrdinary unweighted mass instead obeys $d{\int }_0^{L}qdx/dt=h_1-g_1+\beta [q(L,t)-q(0,t)]$. Replacing the weighted measure by $dy$ or deleting the constant mode would therefore give an incorrect solution.

There is a subtle point when checking nonzero derivative data: every homogeneous eigenfunction has zero endpoint derivative, but differentiating the boundary-forcing sum and its time integral term by term at an endpoint is not valid. The kernel is singular as $t-s\downarrow 0$, and the prescribed derivatives are interior limits of the full solution. An explicit Laplace transform form verifies those limits without this interchange and gives an alternative purely contour-integral version.

For $\mathrm{Re}p>0$, set $\kappa =\sqrt{p+b^2}$ using the principal square root, and defineBoth solve $(p-\mathcal{A})f=0$, with $u_p(0)=v_p(L)=1$, $u_p^{\prime }(0)=v_p^{\prime }(L)=0$. Their weighted Wronskian is the nonzero constantHence the resolvent kernel for Neumann advection-diffusion on an interval isThe first derivative in $x$ has jump $-e^{-\beta y}$, giving $(p-{\mathcal{A}}_x){\mathcal{R}}_p=e^{-\beta y}{\delta }_y$ in the sense of distributional derivatives. If $G(p)={\int }_0^T{e}^{-ps}{g}_1(s)ds$ and $H(p)={\int }_0^T{e}^{-ps}{h}_1(s)ds$, an equivalent transformed solution isIn fact $v_p^{\prime }(0)=-D_p$ and $u_p^{\prime }(L)=e^{-\beta L}{D}_p$, so $Q_x(0,p)=G(p)$ and $Q_x(L,p)=H(p)$ exactly. The inverse Bromwich contour integral of $e^{pt}Q(x,p)/(2\pi i)$ is an alternative final answer. Its poles at $p=0$ and $p=-{\lambda }_n$ reproduce the constant and decaying modes of $K_{\beta }$, respectively; apparent square-root singularities in the normalized hyperbolic factors are removable. This also independently checks the kernel normalization and boundary signs.

As in question 1, the printed dot on $q_0$ denotes its spatial derivative: the two corner conditions are $q_0^{\prime }(0)=g_1(0)$ and $q_0^{\prime }(L)=h_1(0)$. The initial eigenfunction expansion converges to $q_0$ in the weighted space and, for the stated smooth compatible data, has the appropriate classical initial and boundary limits. Differences of solutions with zero data satisfy $\frac{1}{2}d\Vert q{\Vert }_w^2/dt=-{\int }_0^{L}w|q_x{|}^{2}dx$, proving uniqueness in the usual regular solution class. In the limit $\beta \to 0$, $N_0\to L$, ${\varphi }_n\to \cos (n\pi x/L)$, and the ordinary Neumann heat kernel on an interval and its boundary signs are recovered.

Let $T:H^1(U)\to L^2(\partial U)$ be the trace map from the Sobolev trace theorem. The weak Robin problem for a uniformly elliptic operator uses all of $H^1(U)$ as the test space, since a Robin boundary condition does not require a zero trace. DefineThe weak formulation isFor a classical solution, integration by parts produces the boundary term $-{\int }_{\partial U}(\sum a^{ij}{u}_{x_i}{\nu }_j)v$. Substituting the Robin boundary condition replaces it by ${\int }_{\partial U}\beta uv$, giving the displayed identity.