=Paper= {{Paper |id=Vol-2522/paper7 |storemode=property |title=Determination of the Temperature Field in the Depth of the Body Measurement With a Thermal Imager on the Surface |pdfUrl=https://ceur-ws.org/Vol-2522/paper7.pdf |volume=Vol-2522 |authors=Eugenia Homenok,Ivan Ignatiev,Yana F. Ivanova,Dmitriy A. Tarkhov,Andrey D. Yukhnev,Konstantin Zabello }} ==Determination of the Temperature Field in the Depth of the Body Measurement With a Thermal Imager on the Surface== https://ceur-ws.org/Vol-2522/paper7.pdf
                                                  68


    Determination of the Temperature Field in the Depth of
    the Body Measurement With a Thermal Imager on the
                           Surface*

    Eugenia Homenok, Ivan Ignatiev, Yana F. Ivanova, Dmitriy A. Tarkhov[0000-0002-9431-
                  8241]
                        , Andrey D. Yukhnev and Konstantin Zabello

                       Peter the Great St. Petersburg Polytechnic University
                                    dtarkhov@gmail.com



         Abstract. The use of ultrasound local heating of human body tissues for the treat-
         ment of cancer and other diseases is a promising direction in medicine because
         of its minimally invasive and cheap. The application of this method requires pre-
         cise control of the localization of thermal effects in a given area. For this purpose,
         two main methods are used – ultrasonic thermometry and restoration of the ther-
         mal field inside the body by measurements on its surface. In this article, we pre-
         sent the results of our theoretical and experimental researches of the s method.
         The field was reconstructed at depth. This restoration was based on the results of
         measurements of the thermal imager from the surface of the body during high-
         intensity focused ultrasound (HIFU) heating. At the same time, we were forced
         to digitize the graphic data of temperature distribution, which was given by the
         thermal imager. The results of the restoration were compared with the results of
         measurements by sensors which was located at different depths. As a result of
         the comparison, the error was less than 1°C. Subsequently, we took these results
         as exemplary for testing the accuracy of the ultrasonic thermometry (UST) pro-
         gram. The average time difference between the temperature increase obtained by
         the UST program and the recovered data is 1.5 °C. We consider that these results
         are acceptable for further research on biological tissues.

         Keywords: ultrasonic thermometry method, focused ultrasound therapy, recon-
         struction of the temperature field.


1        Introduction



To control the process of thermal effects in the application of focused ultrasound ther-
apy [1] is a usually used method of ultrasonic thermometry [2-6]. To assess the accu-
racy of the ultrasonic testing method, we conducted an experiment in which HIFU heat-
ing of the test object was performed using the thermal imager measured the temperature

*    Copyright 2019 for this paper by its authors. Use permitted under Creative Commons License
     Attribution 4.0 International (CC BY 4.0).
                                            69


on its surface. According to the results of measurements of the temperature field on the
surface of the test object, the field was restored at depth. Results recovery of tempera-
ture by depth were taken the sample of the test program UST. In this paper we present
the results of the temperature field recovery in the test object. We compare these results
with the results of measurements using sensors.


2      Material and methods

The experiment was carried out heating of the tissue-mimicking material (TMM) using
the HIFU transducer. The scanning plane of the ultrasonic transducer is perpendicular
to the axis of the HIFU transducer. The thermal imager was installed at a height of 15
cm from the surface of the test object and recorded a video with thermal images of the
surface of the material during heating and cooling, as well as with the temperature dis-
tribution along a 28mm long line passing through the heating center. To control the
temperature inside the test object, we installed sensors at depths of 2mm, 4mm, and
6mm and at a distance of 3mm from the center of focus.
    To process data from the thermal imager in Matlab, a program was written that al-
lows you to create a graph of the temperature distribution along the line drawn on the
device screen. Work with it is as follows: first, the program is loaded with video data
obtained by the thermal imager and is divided into frames. Next, the coordinates of the
frame area from which the data should be read are set (in this case, the row with the
minimum and maximum values of temperatures and the temperature distribution graph
were set). Then the characters from the string corresponding to the data type specified
in the program are formed into an array. Thus, the output is a table with the maximum
and minimum temperature values from each frame of interest to us. At the same time,
the original graph is cleared of noise so that the program image consists of a white
"background" and a black line. The resulting images contain the temperature and coor-
dinate of each pixel belonging to the graph, that is, an array A[T, i] and length(n) is
formed, where i is the frame number, n is the number of pixels on the x-axis. Since
every 10th frame was taken in this experiment, the time interval between frames is 1 s,
which means that it is possible to plot the temperature change of a fixed x over time
(step – 1 s).
    The problem of restoring the temperature field induced by focused ultrasound from
measurements on the surface of the test object can be solved in the first approximation
in the entire space or in the half-space, neglecting the boundary effects. In this case, the
error is minimal when the source is located in the depth of the material and small time
intervals.
    The first substitution. In this production, the solution is sought throughout the space.
This result can be used to restore the solution in the depth of the test object. According
to [1], the temperature field in the tissues satisfies the equation
                              ∂𝑇
                                   = 𝐾Δ𝑇 − 𝑏𝑇 + 𝑆(𝑡, 𝐱)                                 (1)
                              ∂𝑡

   Here 𝑇 - the thermal field, 𝐾 - the coefficient of thermal diffusivity, 𝑏 - the coeffi-
cient reflecting the heat transfer of blood, 𝑆(𝑡, 𝒙) - ultrasonic heat source.
                                                     70


   According to [1] in cylindrical coordinate system 𝑆(𝑡, 𝒙) = 𝑄(𝑡) 𝑒𝑥𝑝[ − 𝑟 2 /
𝛽𝑟 ] 𝑒𝑥𝑝[ − 𝑧 2 /𝛽𝑧 ]. The result is a solution of the form

T(t, r, z)
       t
             exp[ − b(t − τ)] 𝑒𝑥𝑝[ − r 2 /(4K(t − τ) + βr )] exp[ − z 2 /(4K(t − τ) + βz )]
= c ∫ Q(τ)                                                                                               𝑑𝜏
      0                                   √4K(t − τ) + βz (4K(t − τ) + βr )

   In the case of a short-term ultrasonic pulse at the initial time, the density of the ul-
trasonic source is expressed in terms of the Delta function 𝑄(𝑡) = 𝑞𝛿(𝑡)
                           𝑒𝑥𝑝[ − 𝑏𝑡] 𝑒𝑥𝑝[ − 𝑟 2 /(4𝐾𝑡 + 𝛽𝑟 )] 𝑒𝑥𝑝[ − 𝑧 2 /(4𝐾𝑡 + 𝛽𝑧 )]
      𝑇(𝑡, 𝑟, 𝑧) = 𝑐𝑞
                                                 √4𝐾𝑡 + 𝛽𝑧 (4𝐾𝑡 + 𝛽𝑟 )

   In the first approximation, we can ignore the boundary effects arising from the fact
that the problem is solved in half-space, not in space. In this case, the error is minimal
when the source is located in the depth of the material and small times.
   If a sequence of actions is performed, then, due to the linearity of equation (1), the
total thermal field is represented as the sum of the thermal fields of the sources.
   S substitution. We will solve equation (1) in half-space. Replacement 𝑇(𝑡, 𝒙) =
𝑢(𝑡, 𝒙) 𝑒𝑥𝑝[ − 𝑏𝑡] equation (1) is given as
                                     ∂𝑢
                                           = 𝐾Δ𝑢 + 𝑃(𝑡, 𝐱).                                        (2)
                                     ∂𝑡

Here 𝑃(𝑡, 𝒙) = 𝑆(𝑡, 𝒙) 𝑒𝑥𝑝[ 𝑏𝑡].
   We will solve the problem (2) in the half-space −∞ ≤ 𝑥 < +∞, −∞ < 𝑦 <
+∞, 0 ≤ 𝑧 < +∞. Due to the fact that the thermal conductivity of the air is 25-30 times
                                                                  𝜕𝑢
less than that of the test object, we take as a boundary condition = 0 by 𝑧 = 0. Then,
                                                                  𝜕𝑧
according to [2]
                  𝑡    +∞    +∞    +∞
𝑢(𝑡, 𝑥, 𝑦, 𝑧) = ∫0 ∫−∞ ∫−∞ ∫0           𝑃(𝜏, 𝜉, 𝜂, 𝜁) 𝐺(𝑡 − 𝜏, 𝑥, 𝑦, 𝑧, 𝜉, 𝜂, 𝜁)𝑑𝜏𝑑𝜉𝑑𝜂𝑑𝜁,          (3)
                                          (𝑧−𝜁)2         (𝑧+𝜁)2
                                   𝑒𝑥𝑝[−         ]+𝑒𝑥𝑝[−        ]            (𝑥−𝜉)2 +(𝑦−𝜂)2
where 𝐺(𝑡, 𝑥, 𝑦, 𝑧, 𝜉, 𝜂, 𝜁) =              4𝐾𝑡            4𝐾𝑡
                                                                    𝑒𝑥𝑝[ −                    ].
                                             8(𝜋𝐾𝑡)3/2                            4𝐾𝑡

According to [1] 𝑃(𝑡, 𝑥, 𝑦, 𝑧) = 𝑒𝑥𝑝[ 𝑏𝑡]𝑄(𝑡) 𝑒𝑥𝑝[ − (𝑥 2 + 𝑦 2 )/𝛽𝑟 ] 𝑒𝑥𝑝[ − (𝑧 −
𝑧0 )2 /𝛽𝑧 ]. Here 𝑧0 - depth of focus of ultrasonic action.
                                                               1      𝑡 𝑒𝑥𝑝[𝑏𝜏]𝑄(𝜏) +∞
   Substituting       in    (3),   receive       𝑢(𝑡, 𝑥, 𝑦, 𝑧) =    ∫              ∫−∞ 𝑒𝑥𝑝[ −
                                                          8(𝜋𝐾)3/2 0 (𝑡−𝜏)3/2
 (𝑥−𝜉)2                         +∞        (𝑦−𝜂) 2                         +∞
        ] 𝑒𝑥𝑝[ − 𝜉 2 /𝛽𝑟 ]𝑑𝜉 ∫−∞ 𝑒𝑥𝑝[ −           ] 𝑒𝑥𝑝[ − 𝜂 2 /𝛽𝑟 ]𝑑𝜂 ∫0 (𝑒𝑥𝑝[ −
4𝐾(𝑡−𝜏)                                   4𝐾(𝑡−𝜏)
 (𝑧−𝜁)2               (𝑧+𝜁)2
        ] + 𝑒𝑥𝑝[ −           ]) 𝑒𝑥𝑝[ − (𝜁 − 𝑧0 )2 /𝛽𝑧 ]𝑑𝜁𝑑𝜏
4𝐾(𝑡−𝜏)              4𝐾(𝑡−𝜏)
   The last three integrals are calculated analytically.
                                                                      71


 +∞               (𝑥−𝜉)2                                                   𝜋𝛽 𝐾(𝑡−𝜏)                           𝑥2        +∞
∫−∞ 𝑒𝑥𝑝[ − 4𝐾(𝑡−𝜏)] 𝑒𝑥𝑝[ − 𝜉 2 /𝛽𝑟 ]𝑑𝜉 = 2√𝛽 +4𝐾(𝑡−𝜏)
                                              𝑟
                                                      𝑒𝑥𝑝[ −           ] ∫ 𝑒𝑥𝑝[ −
                                                             𝛽 +4𝐾(𝑡−𝜏) −∞  𝑟                              𝑟

(𝑦−𝜂)2                                                 𝜋𝛽𝑟 𝐾(𝑡−𝜏)                        𝑦2            +∞               (𝑧−𝜁)2
          ] 𝑒𝑥𝑝[ − 𝜂 2 /𝛽𝑟 ]𝑑𝜂 = 2√                                   𝑒𝑥𝑝[ −                        ] ∫0       𝑒𝑥𝑝[ −             ]+
4𝐾(𝑡−𝜏)                                                𝛽𝑟 +4𝐾(𝑡−𝜏)                  𝛽𝑟 +4𝐾(𝑡−𝜏)                         4𝐾(𝑡−𝜏)
             (𝑧+𝜁)2
𝑒𝑥𝑝[ −                 ] 𝑒𝑥𝑝[ − (𝜁 − 𝑧0 )2 /
             4𝐾(𝑡−𝜏)
                𝜋𝐾(𝑡−𝜏)𝛽𝑧                       −(𝑧−𝑧0 )2                          𝑧𝛽𝑧 +4𝐾(𝑡−𝜏)𝑧0
𝛽𝑧 ]𝑑𝜁 = √                       (𝑒𝑥𝑝[                          ] 𝑒𝑟𝑓𝑐[                                        ] +
                𝛽𝑧 +4𝐾(𝑡−𝜏)                    𝛽𝑧 +4𝐾(𝑡−𝜏)                 2√𝛽𝑧 𝐾(𝑡−𝜏)√𝛽𝑧 +4𝐾(𝑡−𝜏)
         −(𝑧+𝑧0 )2                            −𝑧𝛽𝑧 +4𝐾(𝑡−𝜏)𝑧0
𝑒𝑥𝑝[                   ] 𝑒𝑟𝑓𝑐[                                             ])
       𝛽𝑧 +4𝐾(𝑡−𝜏)                   2√𝛽𝑧 𝐾(𝑡−𝜏)√𝛽𝑧 +4𝐾(𝑡−𝜏)

                                 2 +∞ −𝑡 2
   Here 𝑒𝑟𝑓𝑐[ 𝑥] =                ∫ 𝑒 𝑑𝑡 - additional probability integral.
                                √𝜋 𝑥
   As a result, we get:
                       𝛽𝑟 √𝛽𝑧    𝑡               𝑒𝑥𝑝[𝑏𝜏]𝑄(𝜏)
𝑢(𝑡, 𝑥, 𝑦, 𝑧) =                 ∫0 𝛽 +4𝐾(𝑡−𝜏)√𝛽 +4𝐾(𝑡−𝜏) 𝑒𝑥𝑝[ −
                         2               𝑟                 𝑧
  𝑥 2 +𝑦 2                   −(𝑧−𝑧0 )2                           𝑧𝛽𝑧 +4𝐾(𝑡−𝜏)𝑧0
               ](𝑒𝑥𝑝[                          ] 𝑒𝑟𝑓𝑐[                                        ]+
𝛽𝑟 +4𝐾(𝑡−𝜏)        𝛽𝑧 +4𝐾(𝑡−𝜏)        2√𝛽𝑧 𝐾(𝑡−𝜏)√𝛽𝑧 +4𝐾(𝑡−𝜏)
         −(𝑧+𝑧0 )2             −𝑧𝛽𝑧 +4𝐾(𝑡−𝜏)𝑧0
+ 𝑒𝑥𝑝[                    ] 𝑒𝑟𝑓𝑐[                                               ])𝑑𝜏                                               (4)
          𝛽𝑧 +4𝐾(𝑡−𝜏)                        2√𝛽𝑧 𝐾(𝑡−𝜏)√𝛽𝑧 +4𝐾(𝑡−𝜏)

   The surface temperature distribution that can be measured is as follows
                                     𝑡             𝑒𝑥𝑝[𝑏𝜏]𝑄(𝜏)
𝑢(𝑡, 𝑥, 𝑦, 0) = 𝛽𝑟 √𝛽𝑧 ∫0                                                        𝑒𝑥𝑝[ −
                                         (𝛽𝑟 +4𝐾(𝑡−𝜏))√𝛽𝑧 +4𝐾(𝑡−𝜏)
  𝑥 2 +𝑦 2                    −𝑧0 2                            2√𝐾(𝑡−𝜏)𝑧0
               ] 𝑒𝑥𝑝[                         ] 𝑒𝑟𝑓𝑐[                             ]𝑑𝜏                                              (5)
𝛽𝑟 +4𝐾(𝑡−𝜏)              𝛽𝑧 +4𝐾(𝑡−𝜏)                     √𝛽𝑧 √𝛽𝑧 +4𝐾(𝑡−𝜏)

    By measurements using the least-squares method, you can restore the coefficients
𝛽𝑟 , 𝛽𝑧 , 𝐾, 𝑏 and, if necessary, 𝑧0 .
                          𝑞 if 𝜏 ≤ 𝑇
    If 𝑏 = 0, 𝑄(𝜏) = {                :
                          0 if 𝜏 > 𝑇
                                                                  𝑥2 +𝑦2              −𝑧0 2            2√𝐾(𝑡−𝜏)𝑧0
                                                       𝑒𝑥𝑝[−               ] 𝑒𝑥𝑝[             ] 𝑒𝑟𝑓𝑐[                 ]
                           𝑚𝑖𝑛(𝑡,𝑇)                            𝛽𝑟 +4𝐾(𝑡−𝜏)        𝛽𝑧 +4𝐾(𝑡−𝜏)        √𝛽𝑧 √𝛽𝑧 +4𝐾(𝑡−𝜏)
𝑢(𝑡, 𝑥, 𝑦, 0) = 𝛽𝑟 √𝛽𝑧 𝑞 ∫0                                                                                                   𝑑𝜏
                                                                           (𝛽𝑟 +4𝐾(𝑡−𝜏))√𝛽𝑧 +4𝐾(𝑡−𝜏)

   Due to the fact that the measurements are carried out at 𝑡 > 𝑇, receive
                                                        𝑥2 +𝑦2             −𝑧 2
                                                                   0 ] 𝑒𝑟𝑓𝑐[                2√𝐾𝜏𝑧0
                                               𝑒𝑥𝑝[−𝛽 +4𝐾𝜏] 𝑒𝑥𝑝[𝛽 +4𝐾𝜏                                ]
                           𝑡                             𝑟                 𝑧              √𝛽𝑧 √𝛽𝑧 +4𝐾𝜏
𝑢(𝑡, 𝑥, 𝑦, 0) = 𝛽𝑟 √𝛽𝑧 𝑞 ∫𝑡−𝑇                                                                                  𝑑𝜏                  (6)
                                                                 (𝛽𝑟 +4𝐾𝜏)√𝛽𝑧 +4𝐾𝜏


   Introduce the notation 𝐴 = 𝛽𝑟 √𝛽𝑧 𝑞, 𝐵 = 4𝐾, then (6) is slightly simplified
                                              𝑥2 +𝑦2
                                                   0 ] 𝑒𝑟𝑓𝑐[   −𝑧 2               √𝐵𝜏𝑧0
                                 𝑒𝑥𝑝[−𝛽 +𝐵𝜏] 𝑒𝑥𝑝[𝛽 +𝐵𝜏                                     ]
                    𝑡                          𝑟               𝑧                √𝛽𝑧 √𝛽𝑧 +𝐵𝜏
𝑢(𝑡, 𝑥, 𝑦, 0) = 𝐴 ∫𝑡−𝑇                                                                         𝑑𝜏                                  (7)
                                                       (𝛽𝑟 +𝐵𝜏)√𝛽𝑧 +𝐵𝜏

   Let the measurements take place over a period of time 𝛥𝑡 = ℎ, 𝑡 = 𝑛ℎ, 𝑇 = 𝑁ℎ, 𝜏 =
𝑖ℎ. To estimate the integral, we use the formula of averages
                                                  72

                                    𝑥2 +𝑦2            −𝑧0 2            √𝐵𝜏𝑛+𝑖 𝑧0
                           𝑒𝑥𝑝[−             ] 𝑒𝑥𝑝[          ] 𝑒𝑟𝑓𝑐[                ]
                                   𝛽𝑟 +𝐵𝜏𝑛+𝑖       𝛽𝑧 +𝐵𝜏𝑛+𝑖         √𝛽𝑧 √𝛽𝑧 +𝐵𝜏𝑛+𝑖
𝑢(𝑛ℎ, 𝑥, 𝑦, 0) = 𝐴ℎ ∑𝑁
                     𝑖=1                                                            ,    (8)
                                            (𝛽𝑟 +𝐵𝜏𝑛+𝑖 )√𝛽𝑧 +𝐵𝜏𝑛+𝑖

   where 𝜏𝑘 = (𝑘 − 𝑁 − 0.5)ℎ
   Parameters 𝛽𝑟 , 𝛽𝑧 , 𝐴, 𝐵 are found by the least-squares method by minimizing the ex-
pression
∑𝑛,𝑗,𝑘(𝑢(𝑛ℎ, 𝑥𝑗 , 𝑦𝑘 , 0) − 𝑢𝑛 𝑗 )2                                                      (9)
                               𝑘

    where 𝑢𝑛 𝑗 - the measured temperature at the corresponding point.
              𝑘



3       Calculation

Due to the fact that the measurements of the thermal field sensors were carried out in
the test object at a shallow depth, the s substitution of the problem was used. Figure 1
shows the change in the temperature increment in time on the surface of the test object
in the heating center according to the thermal imager and reduced by the formula (6).




Fig. 1. Comparison of the measured thermal imager and the reduced temperature in the heating
                                          center

The graph shows that the difference between the temperature measured by the thermal
imager and the temperature restored by the formula (6) in the heating center by the
above method is about a degree.
                                            73




Fig. 2. Comparison of the measured thermal imager and the restored temperature at a distance
                                of 1.2 mm from the center

It can be seen from figure 2 that the difference between the temperature measured by
the thermal imager and the temperature restored by the formula (6) at a distance of 1.2
mm from the heating center by the above method is about two degrees.




Fig. 3. Comparison of the measured thermal imager and the restored temperature at a distance
                                of 2.4 mm from the center
                                            74


It can be seen from figure 3 that the difference between the temperature measured by
the thermal imager and the temperature restored by the formula (6) at a distance of 2.4
mm from the heating center by the above method is about three degrees.
    Figures 1-3 show that the proposed model of temperature distribution during HIFU
heating of the test object reflects the experimental data quite well. The error increase
while rising the distance from the heating center is caused by the neglect of heat transfer
into the air.
    On fig. 4-6 the change of temperature increment in time at a depth of 2 mm, 4 mm
and 6 mm from the surface of the test object and a distance of 3mm from the heating
center, restored according to the thermal imager according to the formula (6) and meas-
ured by a temperature sensor is given.




Fig. 4. Comparison of the measured and reduced temperature at a depth of 2mm from the sur-
              face of the TMM and a distance of 3m from the heating center

It can be seen from the graph that the difference between the measured temperature and
the thermal imager restored according to the above method is about a degree.
                                           75




 Fig. 5. Comparison of measured and restored temperature at a depth of 4 mm from the TMM
                  surface and a distance of 3 mm from the heating center

It can be seen from the figure that the difference between the measured temperature and
the temperature which was restored by the above method at a depth of 4 mm from the
surface is about 0.5 degrees.




Fig. 6. Comparison of the measured and reduced temperature at a depth of 6m from the TMM
                   surface and a distance of 3m from the heating center

It can be seen from the graph that the difference between the measured temperature and
the thermal imager restored according to the above method is about a half degree.
                                             76


4      Results and Discussion

Based on the above results, we can conclude that the described method of restoring the
temperature field by the results of measurement on the surface can be considered quite
accurate. It should be especially noted that the error in determining the temperature
field in the depth of the test object (the measurements of which were not used in the
construction of the model) is no more than an error in the approximation of the temper-
ature field on the surface (which was used in the construction of the model). In addition,
from a computational point of view, it is more efficient than the UST method. Besides
the UST method, we tested the neural network approach for the described problem [8-
16]. Due to the fact that it is necessary to use sufficiently large neural networks and
labor-intensive training procedures to obtain acceptable accuracy, the specified neural
network approach is not suitable for monitoring the temperature field during HIFU-
heating as the real-time operation is required.


5      Conclusions

The results of the studies that we have given above, allow us to regard the described
method as quite accurate and effective. We recommend using that method in ultrasound
therapy in the treatment of cancer, malignant tumors, and other diseases.


Acknowledgment.
The study was carried out in the framework of Project No. RFMEFI57818X0263
supported by the Education and Science Ministry of the Russian Federation for 2018-
2019.


       References
 1. Khokhlova V.A., Crum L.A., ter Haar G., Aubry J.-F. (Eds.) High Intensity Focused Ultra-
    sound Therapy. Berlin: Springer 500 p. (2017)
 2. Zhou Y. Noninvasive Thermometry in High-Intensity Focused Ultrasound Ablation. Ultra-
    sound Quarterly. 33 (4). 253-260. (2017)
 3. Hsiao Y-S., Deng C.X. Calibration and Evaluation of Ultrasound Thermography using In-
    frared Imaging. Ultrasound Med Biol. 42(2). 503–517. (2016)
 4. Ebbini E.S., Simon C., Liu D., Real-Time Ultrasound Thermography and Thermometry
    IEEE Signal Process. Mag., 35, 166-174 (2018)
 5. Berkovich A.E., Smirnov E.M., Yukhnev A.D., Gataulin Y.A., Sinitsyna D.E., Tarkhov
    D.A. Development of ultrasound thermometry technique using tissue-mimicking phantom.
    IOP Conf. Series: Journal of Physics: Conf. Series 1044 12023 P.7 (2018). DOI:
    10.1088/1742-6596/1044/1/012023
 6. Anand A., Kaczkowski P. Noninvasive measurement of local thermal diffusivity using
    backscatter ultrasound and focused ultrasound heating. Ultrasound in Med. & Biol
 7. A. D. Polyanin. Handbook of Linear Partial Differential Equations for Engineers and Scien-
    tists. Boca Raton-London: Chapman & Hall/CRC Press, 2002
                                              77


 8. Antonov V., Tarkhov D., Vasilyev A. Unified approach to constructing the neural network
    models of real objects. Part 1 // Mathematical Models and Methods in Applied Sciences,
    2018 Volume 41, Issue 18 Pages: 9244-9251
 9. Lagaris I.E., Likas A., Fotiadis D.I. Artificial Neural Networks for Solving Ordinary and
    Partial Differential Equations// IEEE Transactions on Neural Networks. – 1998. – Vol.9,
    No. 5. – pp. 987-1000. DOI: 10.1109/72.712178
10. Dissanayake M.W.M.G., Phan-Thien N. Neural-network-based approximations for solving
    partial differential equations// Communications in Numerical Methods in Engineering. –
    March 1994. – Volume 10, Issue 3. – pp. 195-201.
11. Fasshauer G. E. Solving differential equations with radial basis functions: multilevel meth-
    ods and smoothing// Adv. in Comp. Math. – 1999. – 11. – pp. 139-159
12. Fornberg B., Larsson E. A Numerical Study of some Radial Basis Function based Solution
    Methods for Elliptic PDEs// Computers and Mathematics with Applications. – 2003. – 46.
    – pp. 891-902
13. Galperin E., Pan Z., Zheng Q. Application of global optimization to implicit solution of
    Partial Differential Equations// Computers & Mathematics with Applications. – Pergamon
    Press Ltd. – 1993. – Vol. 25, No. 10/11. – pp. 119-124.
14. Galperin E., Zheng Q. Solution and control of PDE via global optimization methods// Com-
    puters & Mathematics with Applications. – Pergamon Press Ltd. – 1993. – Vol. 25, No.
    10/11. – pp. 103-118
15. Sharan M., Kansa E.J., Gupta S. Application of the Multiquadric method to the numerical
    solution of elliptic partial differential equations// Applied Mathematics and Computation. –
    1997. – 84. – pp. 275-302
16. T.V. Lazovskaya, D.A. Tarkhov and A.N. Vasilyev Parametric Neural Network Modeling
    in Engineering, Recent Patents on Engineering, Volume 11, Number 1, 2017, pp. 10-15