- Email: [email protected]
Fitting of curves and surfaces based on interaction of physical relations and experimental data
RESEARCH NOTE Fitting of curves and surfaces based on interaction of physical relations and experimental data Wojciech Karmowskiand Janusz Orkisz Department of Civil Engineering, (Received July 1982)
Two different matical
approaches
to the solutions
physics are discussed.
second approach Key words:
Technical University of Krakow, Krakow, Poland
is based on experimental
mathematical
of boundary
First, a theoretical
model,
model
value problems
data and a suitable
curve fitting,
boundary
There are two possible approaches to the solution of boundary-value problems in mathematical physics. The first, a theoretical one, assumes an ideal model and carries its analysis out in definite cases. The problem is defined, either by the global formulation if suitable functional is formed, the minimum of which supplies the required solution, or by the local one if physical equations as well as boundary and initial conditions are given. Boundary-value problems formulated in this way are usually solved then either with the aid of analytical or numerical methods, in the latter case mostly the finite element or finite difference methods are being applied. The solution of these equations, although exact for the theoretical model, only approximately describes the physical relationships. Moreover, this inaccuracy is usually magnified by errors resulting from both the numerical methods applied and round off errors. The second approach is based on experimental data and consists of a suitable approximation scheme. In this case only measurement error occurs. However, there are not too many methods that allow experimental data throughout the physical domain to be registered. Such measurements are usually very tedious even when scanning methods are used. Moreover, in this way we can only obtain those quantities, which are directly available from the experiment. Frequently the experiment does not supply us with data of the kind which is directly required in physical or technical analyses. We are forced, therefore, to use at least some of the relationships which hold for the theoretical model. Another difficulty may arise when the experimental data are available for a limited set of points, too scarce in order to obtain the solution of the problem considered. So far, therefore, the theoretical and experimental’73-5 results are obtained separately and only later are they are compared for verification purposes. We propose a combined approach in the present paper. It consists in a simultaneous consideration of all given information, i.e. both theoretical and experimental. The theoretical model is used, to describe general features of the phenomenon as well as
Ltd.
approximation
The
scheme.
value problems
Description of the method
0307-904X/83/01065-05/%03.00 0 1983 Butterworth & Co. (Publishers)
in mathe-
based on an ideal model.
experimental data, that may be considered as a correction factor for the model employed. To obtain the required field quantities one may also use only part of the theoretical equations and part of the experimental data. Remaining data may be used for verification of the results. In the extreme case one may use, as usually is the case, either the theoretical model or the experimental data alone.
Example: 2D case in linear elasticity Let us consider a plane problem of the linear theory of elasticity given in the local formulation as an example of the above approach. The problem is to find stress and displacement fields. The theoretical model is described by the equilibrium and compatibility equations together with Hooke’s law and boundary conditions. On the other hand the plane photo-elasticity techniques allow one to measure isochromatics, isoclinics and isodynes. Experimentally, it is only possible to solve the plane problem when the three fields are known at each point of the domain. In practice the measurements can only be obtained on isolines (isochromatics, isoclinics, isodynes). Unfortunately, the images of these isolines, especially isoclinics are often blurred, their exact position is uncertain, hence large errors are possible. It seems reasonable, therefore, to use both theoretical relations and available experimental data to simultaneously obtain the required solution of the 2D elasticity problem. In the local formulation it is necessary, however, to choose some of these equations in order to obtain a correctly posed boundary-value problem. Thus we have either the equations of the theoretical model:
Appl.
Math.
Modelling,
1983,
Vol. 7, February
65
Research Note
accompanied equations:
by boundary
conditions,
20x, %x - uyy
01-i?*=n
we may always write:
or the experimental
I=; =
k
ox, = m
(2)
ax
ay
0
au,,
au,, -+
-=o
ax
%x=
ay
m
(3)
In equations (l)-(3) n, k, m denote the isochromatic order, isoclinic parameter and isodyne order, respectively. The mixed problem may be formulated in other ways. It is practically important, however, to use the most credible information, e.g. the equilibrium equation from the theory and the isodynes from the experiment. It is possible to consider simultaneously all experimental data and theoretical relationships in the global formulation by constructing a suitable function which should be minimized. The functional may be chosen, e.g. in the form: [(V_o- p)2 + ~(Atru)~]
I = ;
da
I R
(gn -Q)~~L
+;
+;
gl[u&c)--id(
L
+
$ kzl {[al(k) -
@?I - ic(k>2)
2%,(k) hx(k)
-
1 2
-
uyy(k)
ik(k)
where 0 = stress tensor in 2D case: L? = area of domain, p = density of volume force, L = length of boundary; id(k), it(k), ik(k) = value of isodyne, isochromatic and isoclinic at a point ‘k’; Nd, NC, Nk = number of measurements on the isolines, a, b, c,, c,, ck = arbitrarily assumed weighting factors respectively. It is possible to neglect some of the components involved, e.g. dropping physical equations and the compatibility condition (4 = 0). Assuming moreover that c, = ck = 0 we obtain the formulation equivalent to the equations (3) which assures linearization of the problem. In this case isochromatics and isoclinics may be used for verification of the stress field obtained.
Formulation
of functional
66
Appl.
Math.
Modelling,
L
s2
(5)
in which: C2= area of the physical domain; I’ and 0 = differential operators; I’ = unknown function; p = source; L = contour of domain and its length; 4 = value of data on the boundary; a and b = weighting factors connected with experimental quantities, m = value of experimental quantity; N = number of points of measurement. If the weight ‘b’ is negligible, the minimization of the above functional gives the classical solution based on the theoretical approach alone. In the opposite case when weight b is infinite, one obtains the classical method of interpretation of experimental data. Increasing the weight b from zero, one increases the influence of experiment on final results. Usually b must not be increased too much, however, since experimental data may not be sufficient for a reasonable approximation of results, because the value of the unknown function is not defined at points distant from the points of measurement. Only the first component of the function contains all points of the domain. Thus it may be used for extrapolation of data over the whole domain. Interpretation of experimental data is also possible using a general definition of the first component of the function, with the physical structure of the discussed problem not being taken into account. Such a function is useful in cases when either the physical equation does not exist or is not known, e.g. for drawing maps. The function is then defined by the condition of the minimum of curvature, that corresponds to the maximum smoothness of the surface. The curvature of the surface is not a unique quantity and may be defined in different ways. Thus we suggest a new definition that expresses a local average curvature of the surface. It is defined by the square root of the average of the second directional derivative taken in all directions in the ‘xy’ plane:
(6) 0
Here: the v-ordinate of a (x, y) point above the ‘xy’ plane, r-position vector on the xy plane. In the Cartesian coordinates this equation may be written as:
in the general case
The functional contains three components. The first, (ZJ is responsible for conformity of the unknown quantity to equations of the theoretical model. The second, (Zb) accounts for satisfaction of boundary conditions, while the third (I,) ensures agreement of results with experimental data. The first two components may be assumed in the form of one of known functionals in mechanics. In an arbitrary case of a boundary problem given in the local formulation: rV(P)=pforPE,Q
c (OI’--_)2dL
LJ
+$k$, W-m>”
or finally the mixed ones, e.g. in the form:
aax, aa,, -+-=
(I’V-ZQ2dQ+; SL J
andOI/(P)=qforPEL
1983,
Vol.
7, February
This quantity is invariant with respect to rotation of the z axis and satisfies the requirement that only that plane has zero curvature. This is not fulfilled in other known definitions. In the case considered the functional is of the form:
I=;
K2(V)d~+;k~r(~-I$)z I Cl
where vk is a measured value of v.
(8)
Research Note
To find the unknown function V we discretize a suitable functional and minimize it by equalizing its gradient to zero. The equations obtained this way are generally nonlinear and their solution is troublesome. It is easy, however, to solve them in the linear case, like the one discussed above. The choice of the functional being in the form of the sum of squared deviations of real values from virtual ones is arbitrary. The forms assumed here were intended to yield linear algebraic equations. Generally, however, it is not so simple especially when there is a physical reason for choosing the other measure of fitting. The measurements are often made with errors distributed according to a certain statistic which is usually not a square function. If the deviations from the real physical positions of points are not large, then the scattering function can be developed in an incomplete power series with only the first three terms left which yield the square function. It cannot be done, however, in those cases where measurements are made with large errors and the complete scattering function has to be used. The functional is then nonlinear, and it is difficult to find its minimum. This arises for instance in photo-elasticity where isolines are blurred creating real difficulty in finding the correct positions of fringes.
1D case Though the formulation of the problem presented above deals with the 2D case, it may also be applied to the 1D case, provided that the definition for the curvature is given by the formula: (9) where s is the arc length. In the 1D case there are not many variables (nodes). Corrected (smoothed) positions of points may be found by applying the fully nonlinear approach. The method of local variations has been modified and applied here.2 Since it deals with a small number of points, it is very suitable in computer calculations. We will not follow this method for the case of determination of a plane curve. We have to find a smooth curve which passes near the experimental points. Deviation from the experimental position is penalized by a special scattering function. The points are moved to obtain a minimum of the curvature. The function, whose minimum gives the correct positions of these points must therefore contain two components. The first one is responsible for obtaining the maximum smoothness of the curve,6 the second ensures that deviations between the experimental and approximated points are as small as possible. Thus we may write: I=
s L
(~~-Tp)ti
angle formed by the considered point and its two neighbours. Since the curvature is a dimensional quantity, it was multiplied by the sum P + Q = IPI + IQ I (see Figure 1). Finally, the following dimensionless curvature was assumed:
K2=(PQ+P-Q,f’Qp
(11)
(PQ -P.Qj2 with the constant metric functions:
factor /3 equal to 1. In terms of trigono-
1+ cos24
K2Z
(1 - ~0~24)
=p
cos2 f#J
(12)
2 sin44
where: # is an angle between P and Q. The following relations hold: ~~ -+ 00 at 4 + 0 and K = 0 for @J= rr. Since: 3 + cos2@
dK2 -= d@
-2
sinQ,
(1 - cos2q93
@E(O,rr)
(13)
~~ is a monotonic function of 4. Quantity p may be interpreted as a penalty function. The penalty is equal to the scattering function for the argument being the deviation of the virtual point from the experimental one. The weight T should also be discussed. It is chosen as follows: if measurements are strict, then the penalty equals 1 for the zero deviation and 0 for a nonzero deviation, thus:
where: K~ = curvature of starting position of curve; K,in minimum of curvature. Since the weight T should be as small as possible in order to obtain the most smooth curve, we assume that: T=
Kz--K&in
=
(15)
(IO)
where: I= functional; L = arc; K = curvature of curve; p = component responsible for conformity of solution to experimental curve; T = weighting factor, respectively. The curvature of the curve given by a finite set of experimental points should be defined. It was defined as the curvature of the parabola which passes through the given point and two neighbouring ones. A question arises here, how to choose the axis of the local coordinate system. The extreme conditions (straight line, edge) are satisfied only by a coordinate system in which the axis z is a bisector of the
\ Figure 1
Appl.
Math.
Modelling,
1983,
Vol. 7, February
67
Research Note Finally, the approximation is achieved by the following algorithm: we define vector v, components of which are halfwidths of the penalty function in x and y directions respectively. One starts from the first point of the curve and calculates the value of the functional (10) at several points lying in the direction given the vector v. The point for which the function achieves the minimum becomes the new position of the point of the curve. This procedure is continued through all points until deviations of these points are less than prescribed. The method presented above may be applied in cases of 1D and such 2D problems, where experimental data are mesured along isolines. To avoid direct solution of the 2D nonlinear problem, a preliminary handling of experimental data on the subsequent isolines may be done, as a series of separate 1D problems. Only afterwards, when deviations of such newly corrected positions of measured data from their unknown true values would be small enough, do we solve this problem as a 2D but linear one. The preliminary handling of the data aims at correcting the initial positions of measured points lying on isolines so as to obtain curves which are sufficiently smooth. One assumes that these points are situated not far from their true positions. Thus further surface fitting to the given data would be considered as fully 2D but a linear problem, i.e. a search for the minimum of a functional of the type (5) with the second order approximation of the penalty function. The approach presented above is particularly advantageous for the analysis of photo-electric images. The preliminary handling of data consists in smooth curves measured along individual fringes.
f
.
I
x
*”
I,
*
X.I
(16)
As a result of the approximation all nodal values were found with an accuracy of 1%. The 2D Boussinesque problem, whose analytical solution is well known, was also tested using the same technique. Similar quality of the final results was obtained. Further 2D tests are currently being conducted.
68
Appl.
Math.
Modelling,
1983,
Vol.
7, February
,
, f
XI I
1 I”
I
t
l
., I
0
. IX
.
‘X 1
IX
II
* IX+
.
, .
m
,t XI
‘II II
.
+
/
x
/
. I
Figure 2
.
x
Figure 3
.
‘x’
I
.
I
.
“I II
t
”
x
I
f
x/ I
V(x,y) =xz+ya
l
+
Il.
Tests So far, 1D and 2D tests have been carried out separately. Moreover, instead of solving the real physical problems, a computer simulation of experiments was applied. Automatically generated ideal results were then subjected to a randomization process. Curves in 1D were defined by the set of experimental points. The aim of each test was to retrieve the initial ideal curve, starting from the random one and using the approximation method previously described. Verification of the method resulted in the search for the maximal dispersion that allows one to recover the ideal curve. From numerous tests performed we show in Figures 2-4 a circle with exact X, random + and smoothed ‘/ ’ results. In the subsequent figures the dispersion is doubled. Thus one may observe its influence on the quality of the final results. Drawing of maps, which is a relatively simple case of the 2D problem has been tested. The simulation technique with the randomization process included was applied. The maximal dispersion was assumed to be as much as 3%, in order to deal with a linear problem. A regular mesh was applied. In the absence of ‘experimental’ points nodes are denoted by X (Figure 5). The initial surface was assumed to be a paraboloid defined by the formula:
I
I
I
x I
I
x
I
x
.
I
x
, t
Figure 4
.
t
X,
1
I
.
I
.
I
.*
1
.
Research Note
.
.
.
.
.
.
.
.
.
.
.
X
.
X
.
X
.
X
.
.
.
.
X
.
X
.
X
.
X
.
.
X
.
X
.
X
.
X
.
.
.
.
X
.
X
.
X
.
X
.
.
X
.
X
.
X
.
X
.
.
.
.
X
.
X
.
X
.
X
.
.
X
.
X
.
X
.
X
.
.
.
.
X
.
X
.
X
.
X
.
.
.
.
.
.
.
.
.
.
.
Figure 5
Conclusions An approximation method based on experimental data which is physically justified has been presented for 1D and 2D problems. The basic idea of this method is to consider experimental results and a physical model (equations) describing the investigated phenomenon simultaneously. If such equations do not exist or if they are not known, they are replaced by requirement of maximal smoothness of the considered surface or curve. In this paper we have discussed the following aspects: (i) formulation of the problem with physical function of scattering of data (real statistics); (ii) numerical method of the solution of 1D and 2D problems. This method admits nonlinearity of the problem and does not impose any particular class of functions from which a solution is to be found; (iii) simulation techniques verifying the employed approach; (iv), computer programs and numerical tests. The ideas presented in the paper may also be seen differently, namely as a general formulation of boundary value problems in mechanics (physics, engineering, etc.) when all available information, of both a theoretical and an experimental nature, is used simultaneously. In such a case, ‘ill-posed’ problems may arise. Therefore, in the local formulation, a preliminary choice of input data may be
necessary. The most credible equations and experimental data should be retained; In the global formulation, however, one may use all information or only part of it according to the form of the function assumed. Usually the choice is arbitrary and the task may be completed in more than one way. Various generalized functions of mechanics may be applied here for the theoretical part (la, Zb) of function (5), depending on which equations are more credible. Using, for example, the Hu-Washizu function, where no relationships between the displacements, stresses and strains are assumed a priori, and adding a corresponding experimental term I,; instead of arriving at the Euler equations resulting from the exact theoretical model, one finally obtains all the pertinent equations of continuum mechanics modified to account for the experimental data. Though the results of tests made so far were very encouraging much more should be done for verification of these ideas. Further numerical tests are planned, as well as the solution of real life boundary value problems. The direct aim of these efforts is to evaluate the usefulness of the application of all available information or only a part of it at once; to find an optimal choice of weighting factors; to evaluate the effectiveness of the suggested approach; and to choose the best field for applications. Also, because of the newly introduced function some basic mathematical research should be carried out. Finally, it is worth emphasizing that though the approach presented here was given as a method of surface fitting to measured data in the domain of 2D problems of experimental mechanics (photo-elasticity, holography) a similar approach may be used in various boundary value problems of physics and engineering, for example in the interpretation of any interference images and in drawing maps.
References Barnhill, R. D. ‘Representation and approximation of surfaces’. In ‘Math. software III’ (ed. Price, .I.), Academic Press, New York, 1911 Chernoushko, F. L. ‘Variational problems of control and system mechanics’ (Edited) Nauka Press, Moscow, 1973 (in Russian) DeFloriani, L. and Dottori, G. ‘Surface interpolation methods over rectangular and triangular grids’, 2nd Int Congr.Num. Meth. for Eigng, Paris, 1980 Lawson. C. L. ‘Software for C’ surface interuolation’. In ‘Math. Software III’ (ed. Rice, J.), Academic Press,New York, 1977 Schumaker, L. L. ‘Fitting surfaces to scattered data’. In ‘Approximation theory II’ (ed. Lorentz, G. G. et al.), Academic Press, New York, 1976 Zavyalov. Yu. S. ef al. ‘Methods of spline functions’ (Edited) Nauka Press, Moscow, 1980 (in Russian)
Appl.
Math.
Modelling,
1983,
Vol.
7, February
69
Two different matical
approaches
to the solutions
physics are discussed.
second approach Key words:
Technical University of Krakow, Krakow, Poland
is based on experimental
mathematical
of boundary
First, a theoretical
model,
model
value problems
data and a suitable
curve fitting,
boundary
There are two possible approaches to the solution of boundary-value problems in mathematical physics. The first, a theoretical one, assumes an ideal model and carries its analysis out in definite cases. The problem is defined, either by the global formulation if suitable functional is formed, the minimum of which supplies the required solution, or by the local one if physical equations as well as boundary and initial conditions are given. Boundary-value problems formulated in this way are usually solved then either with the aid of analytical or numerical methods, in the latter case mostly the finite element or finite difference methods are being applied. The solution of these equations, although exact for the theoretical model, only approximately describes the physical relationships. Moreover, this inaccuracy is usually magnified by errors resulting from both the numerical methods applied and round off errors. The second approach is based on experimental data and consists of a suitable approximation scheme. In this case only measurement error occurs. However, there are not too many methods that allow experimental data throughout the physical domain to be registered. Such measurements are usually very tedious even when scanning methods are used. Moreover, in this way we can only obtain those quantities, which are directly available from the experiment. Frequently the experiment does not supply us with data of the kind which is directly required in physical or technical analyses. We are forced, therefore, to use at least some of the relationships which hold for the theoretical model. Another difficulty may arise when the experimental data are available for a limited set of points, too scarce in order to obtain the solution of the problem considered. So far, therefore, the theoretical and experimental’73-5 results are obtained separately and only later are they are compared for verification purposes. We propose a combined approach in the present paper. It consists in a simultaneous consideration of all given information, i.e. both theoretical and experimental. The theoretical model is used, to describe general features of the phenomenon as well as
Ltd.
approximation
The
scheme.
value problems
Description of the method
0307-904X/83/01065-05/%03.00 0 1983 Butterworth & Co. (Publishers)
in mathe-
based on an ideal model.
experimental data, that may be considered as a correction factor for the model employed. To obtain the required field quantities one may also use only part of the theoretical equations and part of the experimental data. Remaining data may be used for verification of the results. In the extreme case one may use, as usually is the case, either the theoretical model or the experimental data alone.
Example: 2D case in linear elasticity Let us consider a plane problem of the linear theory of elasticity given in the local formulation as an example of the above approach. The problem is to find stress and displacement fields. The theoretical model is described by the equilibrium and compatibility equations together with Hooke’s law and boundary conditions. On the other hand the plane photo-elasticity techniques allow one to measure isochromatics, isoclinics and isodynes. Experimentally, it is only possible to solve the plane problem when the three fields are known at each point of the domain. In practice the measurements can only be obtained on isolines (isochromatics, isoclinics, isodynes). Unfortunately, the images of these isolines, especially isoclinics are often blurred, their exact position is uncertain, hence large errors are possible. It seems reasonable, therefore, to use both theoretical relations and available experimental data to simultaneously obtain the required solution of the 2D elasticity problem. In the local formulation it is necessary, however, to choose some of these equations in order to obtain a correctly posed boundary-value problem. Thus we have either the equations of the theoretical model:
Appl.
Math.
Modelling,
1983,
Vol. 7, February
65
Research Note
accompanied equations:
by boundary
conditions,
20x, %x - uyy
01-i?*=n
we may always write:
or the experimental
I=; =
k
ox, = m
(2)
ax
ay
0
au,,
au,, -+
-=o
ax
%x=
ay
m
(3)
In equations (l)-(3) n, k, m denote the isochromatic order, isoclinic parameter and isodyne order, respectively. The mixed problem may be formulated in other ways. It is practically important, however, to use the most credible information, e.g. the equilibrium equation from the theory and the isodynes from the experiment. It is possible to consider simultaneously all experimental data and theoretical relationships in the global formulation by constructing a suitable function which should be minimized. The functional may be chosen, e.g. in the form: [(V_o- p)2 + ~(Atru)~]
I = ;
da
I R
(gn -Q)~~L
+;
+;
gl[u&c)--id(
L
+
$ kzl {[al(k) -
@?I - ic(k>2)
2%,(k) hx(k)
-
1 2
-
uyy(k)
ik(k)
where 0 = stress tensor in 2D case: L? = area of domain, p = density of volume force, L = length of boundary; id(k), it(k), ik(k) = value of isodyne, isochromatic and isoclinic at a point ‘k’; Nd, NC, Nk = number of measurements on the isolines, a, b, c,, c,, ck = arbitrarily assumed weighting factors respectively. It is possible to neglect some of the components involved, e.g. dropping physical equations and the compatibility condition (4 = 0). Assuming moreover that c, = ck = 0 we obtain the formulation equivalent to the equations (3) which assures linearization of the problem. In this case isochromatics and isoclinics may be used for verification of the stress field obtained.
Formulation
of functional
66
Appl.
Math.
Modelling,
L
s2
(5)
in which: C2= area of the physical domain; I’ and 0 = differential operators; I’ = unknown function; p = source; L = contour of domain and its length; 4 = value of data on the boundary; a and b = weighting factors connected with experimental quantities, m = value of experimental quantity; N = number of points of measurement. If the weight ‘b’ is negligible, the minimization of the above functional gives the classical solution based on the theoretical approach alone. In the opposite case when weight b is infinite, one obtains the classical method of interpretation of experimental data. Increasing the weight b from zero, one increases the influence of experiment on final results. Usually b must not be increased too much, however, since experimental data may not be sufficient for a reasonable approximation of results, because the value of the unknown function is not defined at points distant from the points of measurement. Only the first component of the function contains all points of the domain. Thus it may be used for extrapolation of data over the whole domain. Interpretation of experimental data is also possible using a general definition of the first component of the function, with the physical structure of the discussed problem not being taken into account. Such a function is useful in cases when either the physical equation does not exist or is not known, e.g. for drawing maps. The function is then defined by the condition of the minimum of curvature, that corresponds to the maximum smoothness of the surface. The curvature of the surface is not a unique quantity and may be defined in different ways. Thus we suggest a new definition that expresses a local average curvature of the surface. It is defined by the square root of the average of the second directional derivative taken in all directions in the ‘xy’ plane:
(6) 0
Here: the v-ordinate of a (x, y) point above the ‘xy’ plane, r-position vector on the xy plane. In the Cartesian coordinates this equation may be written as:
in the general case
The functional contains three components. The first, (ZJ is responsible for conformity of the unknown quantity to equations of the theoretical model. The second, (Zb) accounts for satisfaction of boundary conditions, while the third (I,) ensures agreement of results with experimental data. The first two components may be assumed in the form of one of known functionals in mechanics. In an arbitrary case of a boundary problem given in the local formulation: rV(P)=pforPE,Q
c (OI’--_)2dL
LJ
+$k$, W-m>”
or finally the mixed ones, e.g. in the form:
aax, aa,, -+-=
(I’V-ZQ2dQ+; SL J
andOI/(P)=qforPEL
1983,
Vol.
7, February
This quantity is invariant with respect to rotation of the z axis and satisfies the requirement that only that plane has zero curvature. This is not fulfilled in other known definitions. In the case considered the functional is of the form:
I=;
K2(V)d~+;k~r(~-I$)z I Cl
where vk is a measured value of v.
(8)
Research Note
To find the unknown function V we discretize a suitable functional and minimize it by equalizing its gradient to zero. The equations obtained this way are generally nonlinear and their solution is troublesome. It is easy, however, to solve them in the linear case, like the one discussed above. The choice of the functional being in the form of the sum of squared deviations of real values from virtual ones is arbitrary. The forms assumed here were intended to yield linear algebraic equations. Generally, however, it is not so simple especially when there is a physical reason for choosing the other measure of fitting. The measurements are often made with errors distributed according to a certain statistic which is usually not a square function. If the deviations from the real physical positions of points are not large, then the scattering function can be developed in an incomplete power series with only the first three terms left which yield the square function. It cannot be done, however, in those cases where measurements are made with large errors and the complete scattering function has to be used. The functional is then nonlinear, and it is difficult to find its minimum. This arises for instance in photo-elasticity where isolines are blurred creating real difficulty in finding the correct positions of fringes.
1D case Though the formulation of the problem presented above deals with the 2D case, it may also be applied to the 1D case, provided that the definition for the curvature is given by the formula: (9) where s is the arc length. In the 1D case there are not many variables (nodes). Corrected (smoothed) positions of points may be found by applying the fully nonlinear approach. The method of local variations has been modified and applied here.2 Since it deals with a small number of points, it is very suitable in computer calculations. We will not follow this method for the case of determination of a plane curve. We have to find a smooth curve which passes near the experimental points. Deviation from the experimental position is penalized by a special scattering function. The points are moved to obtain a minimum of the curvature. The function, whose minimum gives the correct positions of these points must therefore contain two components. The first one is responsible for obtaining the maximum smoothness of the curve,6 the second ensures that deviations between the experimental and approximated points are as small as possible. Thus we may write: I=
s L
(~~-Tp)ti
angle formed by the considered point and its two neighbours. Since the curvature is a dimensional quantity, it was multiplied by the sum P + Q = IPI + IQ I (see Figure 1). Finally, the following dimensionless curvature was assumed:
K2=(PQ+P-Q,f’Qp
(11)
(PQ -P.Qj2 with the constant metric functions:
factor /3 equal to 1. In terms of trigono-
1+ cos24
K2Z
(1 - ~0~24)
=p
cos2 f#J
(12)
2 sin44
where: # is an angle between P and Q. The following relations hold: ~~ -+ 00 at 4 + 0 and K = 0 for @J= rr. Since: 3 + cos2@
dK2 -= d@
-2
sinQ,
(1 - cos2q93
@E(O,rr)
(13)
~~ is a monotonic function of 4. Quantity p may be interpreted as a penalty function. The penalty is equal to the scattering function for the argument being the deviation of the virtual point from the experimental one. The weight T should also be discussed. It is chosen as follows: if measurements are strict, then the penalty equals 1 for the zero deviation and 0 for a nonzero deviation, thus:
where: K~ = curvature of starting position of curve; K,in minimum of curvature. Since the weight T should be as small as possible in order to obtain the most smooth curve, we assume that: T=
Kz--K&in
=
(15)
(IO)
where: I= functional; L = arc; K = curvature of curve; p = component responsible for conformity of solution to experimental curve; T = weighting factor, respectively. The curvature of the curve given by a finite set of experimental points should be defined. It was defined as the curvature of the parabola which passes through the given point and two neighbouring ones. A question arises here, how to choose the axis of the local coordinate system. The extreme conditions (straight line, edge) are satisfied only by a coordinate system in which the axis z is a bisector of the
\ Figure 1
Appl.
Math.
Modelling,
1983,
Vol. 7, February
67
Research Note Finally, the approximation is achieved by the following algorithm: we define vector v, components of which are halfwidths of the penalty function in x and y directions respectively. One starts from the first point of the curve and calculates the value of the functional (10) at several points lying in the direction given the vector v. The point for which the function achieves the minimum becomes the new position of the point of the curve. This procedure is continued through all points until deviations of these points are less than prescribed. The method presented above may be applied in cases of 1D and such 2D problems, where experimental data are mesured along isolines. To avoid direct solution of the 2D nonlinear problem, a preliminary handling of experimental data on the subsequent isolines may be done, as a series of separate 1D problems. Only afterwards, when deviations of such newly corrected positions of measured data from their unknown true values would be small enough, do we solve this problem as a 2D but linear one. The preliminary handling of the data aims at correcting the initial positions of measured points lying on isolines so as to obtain curves which are sufficiently smooth. One assumes that these points are situated not far from their true positions. Thus further surface fitting to the given data would be considered as fully 2D but a linear problem, i.e. a search for the minimum of a functional of the type (5) with the second order approximation of the penalty function. The approach presented above is particularly advantageous for the analysis of photo-electric images. The preliminary handling of data consists in smooth curves measured along individual fringes.
f
.
I
x
*”
I,
*
X.I
(16)
As a result of the approximation all nodal values were found with an accuracy of 1%. The 2D Boussinesque problem, whose analytical solution is well known, was also tested using the same technique. Similar quality of the final results was obtained. Further 2D tests are currently being conducted.
68
Appl.
Math.
Modelling,
1983,
Vol.
7, February
,
, f
XI I
1 I”
I
t
l
., I
0
. IX
.
‘X 1
IX
II
* IX+
.
, .
m
,t XI
‘II II
.
+
/
x
/
. I
Figure 2
.
x
Figure 3
.
‘x’
I
.
I
.
“I II
t
”
x
I
f
x/ I
V(x,y) =xz+ya
l
+
Il.
Tests So far, 1D and 2D tests have been carried out separately. Moreover, instead of solving the real physical problems, a computer simulation of experiments was applied. Automatically generated ideal results were then subjected to a randomization process. Curves in 1D were defined by the set of experimental points. The aim of each test was to retrieve the initial ideal curve, starting from the random one and using the approximation method previously described. Verification of the method resulted in the search for the maximal dispersion that allows one to recover the ideal curve. From numerous tests performed we show in Figures 2-4 a circle with exact X, random + and smoothed ‘/ ’ results. In the subsequent figures the dispersion is doubled. Thus one may observe its influence on the quality of the final results. Drawing of maps, which is a relatively simple case of the 2D problem has been tested. The simulation technique with the randomization process included was applied. The maximal dispersion was assumed to be as much as 3%, in order to deal with a linear problem. A regular mesh was applied. In the absence of ‘experimental’ points nodes are denoted by X (Figure 5). The initial surface was assumed to be a paraboloid defined by the formula:
I
I
I
x I
I
x
I
x
.
I
x
, t
Figure 4
.
t
X,
1
I
.
I
.
I
.*
1
.
Research Note
.
.
.
.
.
.
.
.
.
.
.
X
.
X
.
X
.
X
.
.
.
.
X
.
X
.
X
.
X
.
.
X
.
X
.
X
.
X
.
.
.
.
X
.
X
.
X
.
X
.
.
X
.
X
.
X
.
X
.
.
.
.
X
.
X
.
X
.
X
.
.
X
.
X
.
X
.
X
.
.
.
.
X
.
X
.
X
.
X
.
.
.
.
.
.
.
.
.
.
.
Figure 5
Conclusions An approximation method based on experimental data which is physically justified has been presented for 1D and 2D problems. The basic idea of this method is to consider experimental results and a physical model (equations) describing the investigated phenomenon simultaneously. If such equations do not exist or if they are not known, they are replaced by requirement of maximal smoothness of the considered surface or curve. In this paper we have discussed the following aspects: (i) formulation of the problem with physical function of scattering of data (real statistics); (ii) numerical method of the solution of 1D and 2D problems. This method admits nonlinearity of the problem and does not impose any particular class of functions from which a solution is to be found; (iii) simulation techniques verifying the employed approach; (iv), computer programs and numerical tests. The ideas presented in the paper may also be seen differently, namely as a general formulation of boundary value problems in mechanics (physics, engineering, etc.) when all available information, of both a theoretical and an experimental nature, is used simultaneously. In such a case, ‘ill-posed’ problems may arise. Therefore, in the local formulation, a preliminary choice of input data may be
necessary. The most credible equations and experimental data should be retained; In the global formulation, however, one may use all information or only part of it according to the form of the function assumed. Usually the choice is arbitrary and the task may be completed in more than one way. Various generalized functions of mechanics may be applied here for the theoretical part (la, Zb) of function (5), depending on which equations are more credible. Using, for example, the Hu-Washizu function, where no relationships between the displacements, stresses and strains are assumed a priori, and adding a corresponding experimental term I,; instead of arriving at the Euler equations resulting from the exact theoretical model, one finally obtains all the pertinent equations of continuum mechanics modified to account for the experimental data. Though the results of tests made so far were very encouraging much more should be done for verification of these ideas. Further numerical tests are planned, as well as the solution of real life boundary value problems. The direct aim of these efforts is to evaluate the usefulness of the application of all available information or only a part of it at once; to find an optimal choice of weighting factors; to evaluate the effectiveness of the suggested approach; and to choose the best field for applications. Also, because of the newly introduced function some basic mathematical research should be carried out. Finally, it is worth emphasizing that though the approach presented here was given as a method of surface fitting to measured data in the domain of 2D problems of experimental mechanics (photo-elasticity, holography) a similar approach may be used in various boundary value problems of physics and engineering, for example in the interpretation of any interference images and in drawing maps.
References Barnhill, R. D. ‘Representation and approximation of surfaces’. In ‘Math. software III’ (ed. Price, .I.), Academic Press, New York, 1911 Chernoushko, F. L. ‘Variational problems of control and system mechanics’ (Edited) Nauka Press, Moscow, 1973 (in Russian) DeFloriani, L. and Dottori, G. ‘Surface interpolation methods over rectangular and triangular grids’, 2nd Int Congr.Num. Meth. for Eigng, Paris, 1980 Lawson. C. L. ‘Software for C’ surface interuolation’. In ‘Math. Software III’ (ed. Rice, J.), Academic Press,New York, 1977 Schumaker, L. L. ‘Fitting surfaces to scattered data’. In ‘Approximation theory II’ (ed. Lorentz, G. G. et al.), Academic Press, New York, 1976 Zavyalov. Yu. S. ef al. ‘Methods of spline functions’ (Edited) Nauka Press, Moscow, 1980 (in Russian)
Appl.
Math.
Modelling,
1983,
Vol.
7, February
69