c18-2 (779609)
Текст из файла
794Chapter 18.Integral Equations and Inverse TheoryK · D · f = σfMultiplying by D1/2 , we getD1/2 · K · D1/2 · h = σh(18.1.8)where h = D1/2 · f. Equation (18.1.8) is now in the form of a symmetric eigenvalueproblem.Solution of equations (18.1.7) or (18.1.8) will in general give N eigenvalues,where N is the number of quadrature points used. For square-integrable kernels,these will provide good approximations to the lowest N eigenvalues of the integralequation. Kernels of finite rank (also called degenerate or separable kernels) haveonly a finite number of nonzero eigenvalues (possibly none). You can diagnosethis situation by a cluster of eigenvalues σ that are zero to machine precision.
Thenumber of nonzero eigenvalues will stay constant as you increase N to improvetheir accuracy. Some care is required here: A nondegenerate kernel can have aninfinite number of eigenvalues that have an accumulation point at σ = 0. Youdistinguish the two cases by the behavior of the solution as you increase N . If yoususpect a degenerate kernel, you will usually be able to solve the problem by analytictechniques described in all the textbooks.CITED REFERENCES AND FURTHER READING:Delves, L.M., and Mohamed, J.L. 1985, Computational Methods for Integral Equations (Cambridge, U.K.: Cambridge University Press). [1]Atkinson, K.E.
1976, A Survey of Numerical Methods for the Solution of Fredholm IntegralEquations of the Second Kind (Philadelphia: S.I.A.M.).18.2 Volterra EquationsLet us now turn to Volterra equations, of which our prototype is the Volterraequation of the second kind,ZtK(t, s)f(s) ds + g(t)f(t) =(18.2.1)aMost algorithms for Volterra equations march out from t = a, building up the solutionas they go. In this sense they resemble not only forward substitution (as discussedSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).is symmetric.
However, since the weights wj are not equal for most quadraturee (equation 18.1.5) is not symmetric. The matrix eigenvaluerules, the matrix Kproblem is much easier for symmetric matrices, and so we should restore thesymmetry if possible. Provided the weights are positive (which they are for Gaussianquadrature), we can define the diagonal matrix D = diag(wj ) and its square root,√D1/2 = diag( wj ).
Then equation (18.1.7) becomes79518.2 Volterra Equationsin §18.0), but also initial-value problems for ordinary differential equations. In fact,many algorithms for ODEs have counterparts for Volterra equations.The simplest way to proceed is to solve the equation on a mesh with uniformspacing:h≡i = 0, 1, . . . , N,b−aN(18.2.2)To do so, we must choose a quadrature rule. For a uniform mesh, the simplestscheme is the trapezoidal rule, equation (4.1.11):ZtiaK(ti , s)f(s) ds = h 12 Ki0 f0 +i−1XKij fj + 12 Kii fi (18.2.3)j=1Thus the trapezoidal method for equation (18.2.1) is:f0 = g0(1 − 12 hKii )fi = h 12 Ki0 f0 +i−1XKij fj + gi ,(18.2.4)i = 1, .
. . , Nj=1(For a Volterra equation of the first kind, the leading 1 on the left would be absent,and g would have opposite sign, with corresponding straightforward changes in therest of the discussion.)Equation (18.2.4) is an explicit prescription that gives the solution in O(N 2 )operations. Unlike Fredholm equations, it is not necessary to solve a system of linearequations. Volterra equations thus usually involve less work than the correspondingFredholm equations which, as we have seen, do involve the inversion of, sometimeslarge, linear systems.The efficiency of solving Volterra equations is somewhat counterbalanced bythe fact that systems of these equations occur more frequently in practice.
If weinterpret equation (18.2.1) as a vector equation for the vector of m functions f(t),then the kernel K(t, s) is an m × m matrix. Equation (18.2.4) must now also beunderstood as a vector equation. For each i, we have to solve the m × m set oflinear algebraic equations by Gaussian elimination.The routine voltra below implements this algorithm. You must supply anexternal function that returns the kth function of the vector g(t) at the point t, andanother that returns the (k, l) element of the matrix K(t, s) at (t, s). The routinevoltra then returns the vector f(t) at the regularly spaced points ti .#include "nrutil.h"void voltra(int n, int m, float t0, float h, float *t, float **f,float (*g)(int, float), float (*ak)(int, int, float, float))Solves a set of m linear Volterra equations of the second kind using the extended trapezoidal rule.On input, t0 is the starting point of the integration and n-1 is the number of steps of size h tobe taken.
g(k,t) is a user-supplied external function that returns gk (t), while ak(k,l,t,s)is another user-supplied external function that returns the (k, l) element of the matrix K(t, s).The solution is returned in f[1..m][1..n], with the corresponding abscissas in t[1..n].{Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).ti = a + ih,796Chapter 18.Integral Equations and Inverse Theoryvoid lubksb(float **a, int n, int *indx, float b[]);void ludcmp(float **a, int n, int *indx, float *d);int i,j,k,l,*indx;float d,sum,**a,*b;}For nonlinear Volterra equations, equation (18.2.4) holds with the product Kii fireplaced by Kii (fi ), and similarly for the other two products of K’s and f’s. Thusfor each i we solve a nonlinear equation for fi with a known right-hand side.Newton’s method (§9.4 or §9.6) with an initial guess of fi−1 usually works verywell provided the stepsize is not too big.Higher-order methods for solving Volterra equations are, in our opinion, not asimportant as for Fredholm equations, since Volterra equations are relatively easy tosolve.
However, there is an extensive literature on the subject. Several difficultiesarise. First, any method that achieves higher order by operating on several quadraturepoints simultaneously will need a special method to get started, when values at thefirst few points are not yet known.Second, stable quadrature rules can give rise to unexpected instabilities inintegral equations. For example, suppose we try to replace the trapezoidal rule inthe algorithm above with Simpson’s rule. Simpson’s rule naturally integrates overan interval 2h, so we easily get the function values at the even mesh points.
For theodd mesh points, we could try appending one panel of trapezoidal rule. But to whichend of the integration should we append it? We could do one step of trapezoidal rulefollowed by all Simpson’s rule, or Simpson’s rule with one step of trapezoidal ruleat the end. Surprisingly, the former scheme is unstable, while the latter is fine!A simple approach that can be used with the trapezoidal method given aboveis Richardson extrapolation: Compute the solution with stepsize h and h/2. Then,assuming the error scales with h2 , computefE =4f(h/2) − f(h)3(18.2.5)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).indx=ivector(1,m);a=matrix(1,m,1,m);b=vector(1,m);t[1]=t0;for (k=1;k<=m;k++) f[k][1]=(*g)(k,t[1]);Initialize.for (i=2;i<=n;i++) {Take a step h.t[i]=t[i-1]+h;for (k=1;k<=m;k++) {sum=(*g)(k,t[i]);Accumulate right-hand side of linearfor (l=1;l<=m;l++) {equations in sum.sum += 0.5*h*(*ak)(k,l,t[i],t[1])*f[l][1];for (j=2;j<i;j++)sum += h*(*ak)(k,l,t[i],t[j])*f[l][j];a[k][l]=(k == l)-0.5*h*(*ak)(k,l,t[i],t[i]);Left-hand side goes}in matrix a.b[k]=sum;}ludcmp(a,m,indx,&d);Solve linear equations.lubksb(a,m,indx,b);for (k=1;k<=m;k++) f[k][i]=b[k];}free_vector(b,1,m);free_matrix(a,1,m,1,m);free_ivector(indx,1,m);18.3 Integral Equations with Singular Kernels797CITED REFERENCES AND FURTHER READING:Linz, P.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.















