Study of Some Inverse Problems for the Biharmonic

Numerical methods

Trefftz method (TM) – Modified Trefftz method (MTM) The Trefftz method (TM) was first proposed in 1926 . The method can be classified as a boundary-type solution procedure. The main idea of the TM is to use particular solutions 35 as the admissible functions, which satisfy the partial differential equation (PDE) exactly. The numerical effort is required only to approximate the boundary conditions [79]. For example, the general solution, for the two-dimensional Laplace equation in simplyconnected planar domain can express by the series form [75] u(r, θ) = c0 + X∞ n=1 cnr n cos(nθ) + dnr n sin(nθ) (1.4.1) In the conventional Trefftz method, the numerical solution for the two-dimensional Laplace equation in simply-connected planar domain is expressed by linear summation of the following bases [80, 81, 82]. {1,r n cos(nθ),r n sin(nθ), n = 1, 2…} (1.4.2) Recently, Liu [10, 83, 84] had modified the T-complete functions in (1.4.2) by considering the characteristic length of the computational domain is R1 = max0≤θ≤2π r(θ) to stabilize the numerical scheme, ( 1, r R1 !n cos(nθ), r R1 !n sin(nθ), n = 1, 2…) (1.4.3) Hence, the following admissible functions with finite terms truncation m can be used to obtain the Trefftz method (TM) and the modified Trefftz method (MTM), respectively as follows: c0 + Xm n=1 cnr n cos(nθ) + dnr n sin(nθ) (1.4.4) c0 + Xm n=1 cn r R1 !n cos(nθ) + dn r R1 !n sin(nθ) (1.4.5) Definition 1.4.1 (Collocation Method). The collocation method has a great advantage to apply on different geometric shapes, and the simplicity for computer programming. In order to applying the collocation method for approximately solving the equations (1.4.4) and (1.4.5) we seek an approximate solution from a finite-dimensional subspace by requiring that the equations (1.4.4) and (1.4.5) be satisfied at only a finite number of so-called collocation points. We choose m + 1 points 0 ≤ θ0 < θ1, …, < θm ≤ 2π such that the interpolation at these grid points with respect to the subspace Xn is uniquely solvable. For example define θi are the equidistant collocated points on [0, 2π] given by θi = ih, f or i = 0, …, m, and, h = 2π m + 1 . (1.4.6) In both Eq. (1.4.4) and (1.4.5) respectively, for each one there are 2m + 1 unknown coefficients, can be obtained by imposing the different collocated points (1.4.6) with [r0(θi), θi], for i = 1, …, 2m + 1, then the collocation Trefftz method (CTM) ia given by c0 + Xm n=1 cnr0(θi) n cos(nθi) + dnr0(θi) n sin(nθi) = u(r0(θi), θi) (1.4.7) and the modified collocation Trefftz method (MCTM) ia given by c0 + Xm n=1 cn r0(θi) R1 !n cos(nθi) + dn r0(θi) R1 !n sin(nθi) = u(r0(θi), θi) (1.4.8) We obtain a linear equations systems with dimensions n = 2m + 1 can be solved for example by the conjugate gradient method to obtain the coefficients cn, dn, for n ∈ N.

The conjugate gradient method

The classical conjugate gradient method is restricted to solving linear systems Ax = b with matrix A symmetric and positive definite. But recently, many generalizations have been given for solving linear systems that are not symmetric and positive definite. We consider an approach to handling non-symmetric linear systems is to convert them to symmetric linear systems and to then use the original conjugate gradient method [28]. We consider the normal equation A T Ax = A T y such that The matrix of coefficients B = A T A is symmetric and positive definite, and therefore, the conjugate gradient method is applicable. Let x0 be an initial guess for the solution x ∗ = B −1A T y. Define r0 = A T y − Bx0 and s0 = r0. For k > 0, define xk+1 = xk + αk sk, αk = ||rk||2 hBsk, ski rk+1 = A T y − Bxk+1 sk+1 = rk+1 + βk sk, βk = ||rk+1||2 ||rk||2 (1.4.9) with the norm and the inner product are both in L 2 . Theorem 1.4.1. Let B be a self-adjoint matrix. Assume C = I − B is a symmetric positive definite operator. Let xk be generated by the conjugate gradient iteration (1.4.9). Then xk → x ∗ superlinearly: ||x ∗ − xk|| ≤ (ck) k ||x ∗ − x0||, k ≥ 0 (1.4.10) with ck → 0, as k → ∞. 1.4.3 Numerical integration Numerical integration formula, or quadrature formula, are methods for the approximate evaluation of definite integrals. Definition 1.4.2. In general, a quadrature formula is a numerical method for approximating an integral of the form Q(f) = Z b a f(x)dx (1.4.11) of a continuous function f over the interval [a, b] with a < b by a weighted sum Qn(f) = Xn j=0 a n j f(x n j ) with quadrature points x n 0 , …, x n n ∈ [a, b] and real quadrature weights a n 0 , …, a n n . 1.4.4 Trigonometric interpolation One important method of constructing degenerate kernels that approximate a given continuous kernel is the interpolation, for which we recall the following theorems.  Definition 1.4.3. For n ∈ N we denote by Tn the linear space of trigonometric polynomials q(t) = Xn m=0 an cos mt + Xn m=1 bm sin mt (1.4.12) with real (or complex) coefficients a0, …, an and b0, …, bn. A trigonometric polynomial q ∈ Tn is said to be of degree n if |an| + |bn| > 0. Theorem 1.4.2. Given 2n + 1 distinct points t0, …, t2n ∈ [0, 2π) and 2n + 1 values y0, …, y2n ∈ R, there exists a uniquely determined trigonometric polynomial qn ∈ Tn with the property qn(tj) = yj , j = 0, …, 2n. Proof. [42]  Theorem 1.4.3. Let tj = jπ n , j = 0, …, 2n − 1, be an equidistant subdivision of the interval [0, 2π] with an even number of grid points. Then, given the values y0, …, y2n−1, there exists a unique trigonometric polynomial of the form qn(t) = α0 2 + Xn−1 k=1 [αk cos kt + βk sin kt] + αn 2 cos nt (1.4.13) satisfying the interpolation property qn(tj) = yj , j = 0, …, 2n − 1 . It’s coefficients are given by αk = 1 n 2 Xn−1 k=0 αk cos kt, k = 0, …, n, βk = 1 n 2 Xn−1 k=0 αk sin kt, k = 0, …, n − 1, Proof. See [41, 42]  Remark 1.4.1. As consequence from theorem 1.4.3. The Lagrange basis for the trigonometric interpolation has the form (see [41]) Lj(t) = 1 2n    1 + 2 Xn−1 k=1 cos k(t − tj) + sin n(t − tj)    (1.4.14) for t ∈ [0, 2π] and j = 0, …, 2n − 1. 1.4.5 The Nyström method The Nyström method or quadrature method, was originally introduced to handle approximations based on numerical integration of the integral operator, the resulting solution is found first at the set of quadrature node points, and then it is extended to all points in Ω by means of a special, and generally quite accurate, interpolation formula. The numerical method is much simpler to 38 implement on a computer, but the error analysis is more sophisticated. We choose a convergent sequence Qn of quadrature formula for the integral in definition 1.4.2. With quadrature points x n 0 , …, x n n ∈ [a, b] and real quadrature weights a n 0 , …, a n n . We approximate the integral operator (Tϕ)(x) = Z b a K(x, y)ϕ(y)dy, x ∈ [a, b]; with continuous kernel K by a sequence of numerical integration operators (Tnϕ)(x) = Xn j=0 ajK(x, yj)ϕ(yj), x ∈ [a, b]; i.e., we apply the quadrature formula for g = K(x, .)ϕ. Then the solution to the integral equation of the second kind is approximated by the solution of ϕn − Tnϕn = f which reduces to solving a finite-dimensional linear system. Theorem 1.4.4. Let ϕn be a solution of ϕn(x) − Xn i=1 Kn(x, yi)ϕn(yi) = f(x), s ∈ [a, b]. Then the values ϕ n j = ϕn(xj), j = 1, …, n, at the quadrature points satisfy the linear system ϕ n j − Xn i=1 Kn(x, yi)ϕ n i = f(xj), j = 1, …, n. Theorem 1.4.5. Assume the quadrature formulas (Qn) are convergent. Then the sequence (Tn) is collectively compact and pointwise convergent (i.e., Tnϕ −→ Tϕ, n → ∞, for all ϕ ∈ C([a, b])), but not norm convergent. Corollary 1.4.1. For a uniquely solvable integral equation of the second kind with a continuous kernel and a continuous right-hand side, the Nyström method with a convergent sequence of quadrature formulas is uniformly convergent. Remark 1.4.2 (see [41]). We will now describe the application of Nyström’s method for the approximate solution of integral equations of the second kind with weakly singular kernels of the form (Tϕ)(x) = Z b a H(|x − y|)K(x, y)ϕ(y)dy, x ∈ Ω; H is continuous and satisfies |H(t)| ≤ Mtα−m for all t > 0 and some positive constants M and α. Then we approximate the weakly singular integral operator by a sequence of numerical integration operators (Tnϕ)(x) = Xn j=1 ajK(x, yj)ϕ(yj) We confine ourselves to a special case by considering a weakly singular operator with a logarithmic singularity (Tϕ)(t) = Z 2π 0 ln(4 sin2 t − τ 2 )K(t, τ)ϕ(τ)dτ, 0 ≤ t ≤ 2π; (1.4.15) 39 in the space C(0, 2π) ⊂ C(R) of 2π-periodic continuous functions. The kernel function K is assumed to be continuous and 2π-periodic with respect to both variables. We construct numerical quadratures for the improper integral (Qϕ)(t) = Z 2π 0 ln(4 sin2 t − τ 2 )ϕ(τ)dτ, 0 ≤ t ≤ 2π; (1.4.16) Using the Lagrange basis 1.4.14 we obtain that [41] (Qnϕ)(t) = 2 Xn−1 j=0 ajR n j (t)ϕ(tj), 0 ≤ t ≤ 2π; (1.4.17) with the quadrature weights, for j = 0, …, 2n − 1. R n j (t) = − 1 n    Xn−1 k=1 1 k cos k(t − tj) + 1 2n cos n(t − tj)    , (1.4.18) and 2 Xn−1 j=0 |R n j (t)| ≤ π √ 2 Remark 1.4.3. Because of singularities in equation (1.4.15) at the two points t = 0 and t = π, we discretizing the equations with equidistant points on [0, 2π] would lead to a poor accuracy. For this reason, it is more appropriate to use a mesh that is graded towards the intersection points. Such a grading can be achieved most efficiently by using a sigmoidal transformation, i.e., a strictly monotonically increasing function as the follows [19]. Definition 1.4.4 (Sigmoidal Transformation). [43] We define the quadrature weights given by wp : [0, 2π] → [0, 2π] wp(t) = 2π [υ(t)]p [υ(t)]p + [υ(2π − t)]p , 0 ≤ t ≤ 2π (1.4.19) with υp is a cubic polynomial given by υp(t) = ( 1 p − 1 2 )  π − t π 3 + 1 p ( t − π π ) + 1 2 the function wp is strictly monotonically increasing, with derivatives vanishing up to a certain order p − 1, p ≥ 2, at the two intersection points t = 0 and t = π with w 0 p (π) = 2, υp(0) = 0, υp(2π) = 1. Remark 1.4.4. The parameter p in the substitution functions is the so-called grading parameter. For larger values of p the grid points are more densely accumulated at the end points of the integration interval (see [11]). To this end we take the weights aj , and the mesh tj as : aj = w 0 p (tj), tj = wp(tj), j = 0, …, 2n − 1.