# Heuristic Search For Physics-Based Problems: Angry Birds in PDDL+

Wiktor Piotrowski1, Yoni Sher1, Sachin Grover1, Roni Stern1,2, Shiwali Mohan1

###### Abstract

This paper studies how a domain-independent planner and combinatorial search
can be employed to play Angry Birds, a well established AI challenge problem.
To model the game, we use PDDL+, a planning language for mixed
discrete/continuous domains that supports durative processes and exogenous
events. The paper describes the model and identifies key design decisions that
reduce the problem complexity. In addition, we propose several domain-specific
enhancements including heuristics and a search technique similar to preferred
operators. Together, they alleviate the complexity of combinatorial search. We
evaluate our approach by comparing its performance with dedicated domain-
specific solvers on a range of Angry Birds levels. The results show that our
performance is on par with these domain-specific approaches in most levels,
even without using our domain-specific search enhancements.

## Introduction

Angry Birds, a wildly popular mobile game, is an open challenge problem for AI
(Renz et al. 2019) that requires reasoning about sequential actions in a
continuous world with discrete exogenous events. Different versions of the
game have proven to be NP-Hard, PSPACE-complete, and EXPTIME-hard (Stephenson,
Renz, and Ge 2020), and the reigning world champion is still a human. The
game’s simple layout and mechanics teamed with human players’ innate spatio-
temporal reasoning and forward state prediction makes for a challenging task.
However, where humans excel, AI agents struggle. Small errors in initial
assumptions can result in scores and states drastically different from the
player’s predictions. Angry Birds requires a holistic understanding of each
individual level and its relevant characteristics. Understanding the relevance
of all features and the sum of all of its parts is not a common trait of AI
approaches. Angry Birds is a difficult and fascinating challenge for
autonomous agents, solving it would prove a significant milestone in AI.

In this work, we present the first successful game playing agent for Angry
Birds that uses a domain-independent planner and combinatorial search. Most
existing planning languages, such a STRIPS (Fikes and Nilsson 1971), PDDL
(McDermott et al. 1998), and PDDL2.1 (Fox and Long 2003), lack features such
as exogenous activity to capture the Angry Birds games. In this paper, we
study how the game mechanics can be encoded using PDDL+ (Fox and Long 2006), a
rich planning language that is designed for mixed discrete/continuous domains.
PDDL+ enables modeling the physics-based dynamics such as effect of gravity,
collision between objects, and trajectory of the bird using equations of
motion and their computational approximations (e.g., Bhaskara’s sine
approximation formula).

Our first contribution is a model of Angry Birds using PDDL+, which includes
several important modeling choices that ensure solution in reasonable time. We
employ time discretization techniques to reason about continuous aspects of
the domain which presents the challenge of combinatorial search. The second
contribution of this work is the development of several search enhancements
designed specifically to constrain solution space in the domain. These
enhancements include domain-specific heuristics and a “preferred states”
mechanism similar to the preferred operators technique (Richter and Helmert
2009). Third contribution is through a comprehensive evaluation against some
of the state-of-the-art agents for Angry Birds (Borovicka, Spetlik, and Rymes
2014; Wang 2017). Our evaluation demonstrates that our techniques can solve a
greater diversity of Angry Bird levels compared to other agents thus,
showcasing that proposed domain-independent search strategies and domain-
specific heuristics make search more efficient. Please note, this is a
preliminary exploratory work incorporating the initial modeling and searching
techniques for discrete-continuous domains, such as basic features of birds of
point and shoot. Since, submitting this work we have expanded the capabilities
to incorporate some of the special features, such as, second tap to initialize
special power of birds, that will be part of the forthcoming submissions.

The paper is organized as follows. We begin by providing relevant background
on the Angry Birds game and PDDL+. This is followed by an analysis of
different planning formalisms and their limitations for this domain. Next, we
present the design of our PDDL+ model for Angry Birds. Then, we describe how
this PDDL+ model is used within a game playing agent, and propose several
domain-specific search enhancements. Finally, we close with an experimental
evaluation of our approach on a variety of Science Birds problems comparing
our results with state-of-the-art baseline agents from previous competitions.
For the convenience of the readers we have provided the complete domain, a
problem file and the solution generated as part of the supplementary material
with the submission.

## Background

Figure 1: Trivial Angry Birds level. Figure 2: Difficult Angry Birds level.

PDDL+ (Fox and Long 2006) is an extension of the well-known Planning Domain
Description Language (PDDL) (McDermott et al. 1998), that maps the planning
constructs to a Hybrid Automata (Henzinger 2000) for improved expressiveness
with numeric fluents and metrics, instantaneous actions, exogenous events and
durative changes called processes. PDDL+ is designed to model and plan for
problems requiring both discrete mode switches, continuous flows with
exogenous environmental changes. It builds on the expressiveness of its
predecessors modeling languages, encapsulating the entire set of features from
PDDL2.1 (Fox and Long 2003) and adding timed-initial literals from PDDL2.2
(Edelkamp and Hoffmann 2004).

In PDDL+ the domain $\mathcal{D}$ for the game environment is a tuple $\langle
F,R,A,E,P\rangle$ where $F$ is a set of discrete state variables; $R$ is a set
of numeric state variables; $A$ is a set of actions the agent can execute; $E$
is a set of instantaneous events triggered in the environment directly or
indirectly due to actions of the agent; $P$ is a set of durative process that
may be active in the environment or triggered by the agent. A state of the
environment $s$ is a complete assignment of values to the state variables
$F\cup R$. Please note, that the domain of numeric state variables does not
include $\bot$ value, i.e., none of the variables are undefined in the state
as the environment and the transitions are completely deterministic. A
planning problem $\Pi=\langle\mathcal{D},s_{I},G\rangle$, where $s_{I}$ is the
initial state, and $s_{I}\in S$ the set of possible states in the environment;
$G$ is the goal condition for the environment.

#### The Angry Birds Game

consists of several different levels. Every level has some allocated numbers
of birds that are launched using a slingshot in an attempt to hit pigs sitting
inside structures built using platforms and blocks. Pigs can be killed by a
direct hit from a bird, or indirectly by falling blocks, explosions, or
falling from a height. Every level consists of different arrangement of
structures and a number of pigs to be killed. Platforms are indestructible
floating objects, while blocks can be destroyed or toppled over. Blocks come
in various shapes, sizes, and materials. For example, blocks made of ice break
easily, stone blocks are harder to break, and TNT crates (shown in fig. 2)
explodes on impact, destroying close objects and launching others in the air.
Some birds have special abilities, activated by tapping on the screen during
flight. While these special abilities are needed to pass some levels, we do
not model them in this paper for simplicity.

Aim of the game is to kill the pigs using minimum number of birds and maximize
the game score. A well-known strategy of angry birds is to hit the weakest
point of the structure, to maximize the impact on the pigs and kill them
quickly. The composition of each level can drastically affect the choice of
strategy for playing the game. For example, Figure 1 shows a simple level with
red birds and a single pig protected by a weak structure. To pass this level,
it is sufficient to shoot a bird directly towards the pig, causing the
structure to collapse on it. In contrast, Figure 2 shows a more difficult
level where the only viable strategy is executing precise shots that set off
an explosive chain reaction. Thus, an intelligent agent playing Angry Birds
needs to be able to predict future states of the world that extend beyond the
immediate and direct consequences of the agent’s actions. As the game is
solved in discrete and continuous state space, it has become the game of
choice for a long running yearly competition, AI Birds (Renz et al. 2015),
organized by Australian National University at the IJCAI conference.

## Modeling Angry Birds

To create a domain for Angry Birds, we need to model the flight of the bird,
collisions between structures, explosions, and structure collapse after
collisions and explosions. In this section we look at these aspects in some
detail. First, we begin with a brief discussion on why we chose PDDL+ as our
modeling language, then describe the PDDL+ components used to model different
objects in the game, and finally the PDDL+ components used to model the game
dynamics.

| PDDL+ | PDDL2.1 | PDDL2.2
---|---|---|---
Bird Flight | $\checkmark$ | $\checkmark$ |
Collisions | $\checkmark$ |  | $\checkmark$
Explosions | $\checkmark$ |  | $\checkmark$
Structure collapse | $\checkmark$ |  | $\checkmark$
Table 1: Language support for crucial Angry Birds features. FSTRIPS+ similar to PDDL2.1 is also able to model only the Bird Flight. | Predicate | Domain | Comments
---|---|---|---
Birds | bird_type | $\mathbb{Z}$ | Different types of birds such as 0 is for the Red bird
x_bird, y_bird | $\mathbb{R}$ | XY location of the bird
v_bird, vx_bird, vy_bird | $\mathbb{R}$ | Cumulative and XY component of velocity of of the bird
m_bird | $\mathbb{R}^{+}$ | Mass of the bird.
bird_id | $\mathbb{Z}$ | Sequence in which the birds will be fired.
Pigs | x_pig, y_pig | $\mathbb{R}$ | XY location of the pig
r_pig, m_pig | $\mathbb{R}^{+}$ | Radius & Mass of the pig
Blocks | x_block, y_block | $\mathbb{R}$ | XY location of the block
block_height, block_width | $\mathbb{R}^{+}$ | Height and Width of the block
block_mass | $\mathbb{R}^{+}$ | Mass of the block
Platforms | x_platform, y_platform | $\mathbb{R}$ | XY location of the platform
platform_width | $\mathbb{R}^{+}$ | Width of the platform.
platform_height | $\mathbb{R}^{+}$ | Height of the platform.
Table 2: Modeling different objects in Angry Birds.

Different aspects of Angry Birds can be modelled with different levels of
expressiveness. For example, PDDL2.1 (Fox and Long 2003) can be used to model
the flight of the bird using numeric fluents similar to an FSTRIPS+ (Ramirez
et al. 2017). Table 1 provides a comparison of PDDL2.1, PDDL2.2 (Edelkamp and
Hoffmann 2004), FSTRIPS+, and PDDL+, and what parts of the game can be
modelled using them. PDDL2.2 provides a close competition for modeling,
however, PDDL+ provides support for all parts of the game.

PDDL+ planning problems are difficult to solve due to immense search spaces,
and complex system dynamics. Indeed, planning in hybrid domains is challenging
because, apart from the state space explosion caused by discrete state
variables, the continuous variables cause the reachability problem to become
undecidable (Alur et al. 1995). Thus, clever design choices when modeling in
PDDL+ are often crucial to solving PDDL+ problems, arguably more so than for
any other planning domain definition languages.

### Angry Birds Objects in PDDL+

The objects in the Angry Birds domain are separated into four types: birds,
pigs, blocks, and platforms. Note that we did not list the slingshot, which is
where the birds are launched from, as an object in the domain. Instead of
explicitly modeling the slingshot, we assign the coordinates of the at-rest
slingshot to every bird that is about to be launched. Avoiding modeling
unnecessary objects helps in maintaining a light-weight domain definition.
Table 2 describes different predicates that were used to model birds, pigs,
blocks and platforms. The state is defined using the XY location of each
object and their specific properties, e.g., mass or radius.

### Angry Birds Dynamics in PDDL+

The PDDL+ domain of Angry Birds features a variety of dynamics which dictate
the change and evolution of the system. This, however, creates planning
problems with vast search spaces and large branching factors, making them
computationally difficult to solve. Thus, mitigating the issues of solvability
and efficiency was at the forefront of the domain’s development at each step
of the process, and is evident in the resulting PDDL+ model’s composition.

#### Launching a Bird

Our PDDL+ model of Angry Birds contains a single action – launching a bird
from the slingshot at a chosen angle. Launching of the bird is split into
three separate phases: selecting the launch velocity, finding the launch
angle, and releasing the bird. In our model, we fix the launch velocity to its
maximum possible value to decrease the solution search space. This decision is
motivated by the fact that maximum velocity shots provide widest range of
targets and the final impact velocity of the bird is directly proportional to
damage it causes to the objects. Moreover, the same level of reachability can
be modeled by changing the angle of slingshot alone instead of modulating the
speed and angle together.

Next, we address the construction of the action responsible for selecting the
launch angle. In a naive approach, one could encode a set of instantaneous
actions that increase or decrease the angle, supplemented by an action
releasing the bird from the slingshot. In this approach, the search space will
include many cycles, corresponding to increasing and decreasing the launch
angle by equal amounts. Instead, we pair the release action with a supporting
process that continuously adjusts the angle as soon as a new bird is placed on
the slingshot, ready to be launched. This follows the Theory of Waiting
(McDermott 2003), which sees the agent idly waiting until the world evolves
into a favorable state in which to execute actions. This encoding reduces the
number of decision points, cycles, and branching factor, mitigating state
space explosion.

Concretely, the model includes an _increase_angle_ process for increasing the
launch angle and a _release_bird_ action (fig. 3) whose effects assign values
to the vertical and horizontal velocity variables based on the angle. The
velocities are then used to model the ballistic flight of the birds. Once a
bird is launched, another process is triggered. This flying process (fig. 4)
models ballistic flight of the active bird, updating its velocity and location
over time, according to the governing equations of motion. Note that this
process represents the domain’s non-linear behavior, which is a major obstacle
in PDDL+ planning as many planners struggle to handle complex dynamics.

     (:action release_bird
        :parameters (?b - bird)
        :precondition (and
            (= (active_bird) (bird_id ?b))
            (not (angle_adjusted))
            (not (bird_released ?b)) )
        :effect (and
            (assign (vy_bird ?b)
                (* (v_bird ?b) SIN(angle)))
            (assign (vx_bird ?b)
                (* (v_bird ?b) COS(angle)))
            (bird_released ?b) (angle_adjusted)) )

Figure 3: PDDL+ action for launching a bird. Since trigonometric functions are
not included in PDDL+, model relies on Bhaskara’s small-angle approximations
for SIN and COS.

#### Events

All the objects interact with the environment and other objects through
events. Below, we discuss three different classes of events for modeling
collisions, motion of the blocks, and auxiliary events at the start and end of
each shot. These events model complex behaviors of destruction, explosion and
structure collapse when hit by birds. In our planning model, all the post
launch system dynamics are caused by events interacting with the effects of
the flying process.

     (:process flying
        :parameters (?b - bird)
        :precondition (and
            (bird_released ?b)
            (= (active_bird) (bird_id ?b))
            (> (y_bird ?b) 0) )
        :effect (and
            (decrease (vy_bird ?b)
                (* #t (gravity) ))
            (increase (y_bird ?b)
                (* #t (vy_bird ?b)))
            (increase (x_bird ?b)
                (* #t (vx_bird ?b)))) )

Figure 4: PDDL+ process of a bird’s flight.

Collisions are the first class of events that are highlighted in an Angry
Birds problem. The main aim of the game is to kill pigs by directly colliding
with them, causing nearby explosions, or destroying blocks such that pigs are
killed by collapsing structures or falling from heights. Each of those cases
is a crucial tactic in Angry Birds, and any accurate model of the game must
encompass such mechanics.

As an example event, direct interactions between a bird and a pig are modeled
based on elastic sphere collisions in two dimensions. Pigs are very fragile
and disappear upon almost any contact from another object. Thus, after a
collision, only the bird’s resulting velocity needs to be calculated,
according to the angle-free elastic collision formula:

$\vec{v}_{b}^{\prime}=\vec{v}_{b}-\frac{2m_{p}}{m_{b}+m_{p}}\frac{\langle\vec{v}_{b}-\vec{v}_{p},\vec{x}_{b}-\vec{x}_{p}\rangle}{||\vec{x}_{b}-\vec{x}_{p}||^{2}}(\vec{x}_{b}-\vec{x}_{p})$
(1)

where $\vec{v}_{b}$ and $\vec{v}_{p}$ are the velocity of the bird and pig,
respectively, $\vec{x}_{b}$ and $\vec{x}_{p}$ are their respective positions,
$m_{b}$ and $m_{p}$ are their respective masses, and the angle brackets denote
an inner product of two vectors. Translating this equation into a PDDL+ event
allows us to model behavior beyond the obvious solutions (i.e., aiming
directly at an exposed pig). It enables the planner to find highly innovative
solutions such as deliberately and directly killing two pigs with one bird
(fig. 6). Similarly, we model events for collisions with the grounds, which
also enables finding trick shots in which the bird bounces off the ground to
hit a pig or a block.

The ground in Angry Birds is often overlooked when playing the game but it can
prove useful when aiming at difficult targets. As such, an event (seen in fig
5 is defined in the domain that influences the bird’s trajectory after
bouncing off the ground, based on energy lost during the collision (i.e.,
ground damper). Figure 7 shows how a ground-bounce event helps in killing a
pig when a direct hit is impossible.

     (:event collision_ground
      :parameters (?b - bird)
      :precondition (and
         (= (active_bird) (bird_id ?b))
         (<= (y_bird ?b) 0)  )
      :effect (and
         (assign (y_bird ?b) 1)
         (assign (vy_bird ?b)
         (* (* (vy_bird ?b) -1)(ground_damper)))
         (assign (bounce_count ?b)
         (+ (bounce_count ?b) 1))) )

Figure 5: Event modeling a bird bouncing off the ground.

Figure 6: Two pigs-one bird plan in execution.

Figure 7: Solution requiring bouncing off the ground and collapsing a
structure to kill the target pig.

Most interactions in a level occurs between birds and blocks in the scene.
Blocks form structures which act as a protection barrier for targeted pigs.
Thus, modeling collisions with blocks is vital to playing Angry Birds.
However, similarly to pig collisions, we opt to only define the event’s
effects in terms of the impact it has on the bird’s velocity and the blocks
stability and life values. There are two events, modeling bird-block
collisions with the distinction based on whether the block is stable enough to
withstand the impact without being knocked out of place or destroyed entirely.
In such cases, the bird bounces back off the block, otherwise it continues
moving forward with a diminished velocity.

Modeling the motion of blocks after impact and chains of block-block
collisions would be prohibitively difficult, requiring a process for tracking
each individual block’s change in position and rotation over time, as well as
a set of events to account for secondary interactions. These additions would
only marginally improve the domain’s accuracy, though the planner would
experience a drastic drop in performance, having to keep track of dozens of
simultaneous non-linear processes in every state after a collision. Instead,
we model one central block-block interaction, namely dislodging or destroying
blocks supporting larger structures. In such cases, an event is triggered for
every supported block, reducing its stability value to 0, repositioning the
block to the ground, and adjusting its life value to account for fall damage.
Analogously, we defined another event for killing pigs which sit atop
collapsing block structures.

The two final block-related events concern exploding TNT crates. Any impact
upon the crate causes an explosion destroying all objects in its immediate
vicinity and displacing objects further away. In our PDDL+ model, we
incorporate two simple events modeling the explosive destruction of pigs and
blocks near a TNT crate. As previously indicated, we omit movement of objects
caused by explosions.

Platforms are indestructible objects which the agent needs to reason with to
avoid poor results. However, almost no meaningful interactions occur between
birds and platforms, and thus our domain contains an event modeling this type
of collision but only to discourage the agent from shooting at platforms. The
event deters the agent by nullifying the bird’s velocity and expiring it
without any impact on the level.

Finally, we encode a supporting event which models the end of life for an
active bird. In the game, a bird expires when it moves beyond the limits of
the scene or when it stops moving. Modeling the exact timing of the end of a
bird’s life is quite difficult. Instead, we define the birds usefulness in
terms of the number of bounces off other objects. Each bounce partially
reduces the bird’s momentum, and its collision impact becomes negligible after
three bounces. Therefore, we define a termination event for the active bird
once it has bounced three times. The result of this event is that the next
bird is placed on the slingshot for launch.

Goals While the goal of Angry Birds is to destroy all the pigs, we consider a
simpler formulation for individual planning problems. The simplified problem
splits the original scenario into single-bird episodes with the goal of
killing at least one pig. During evaluations we start with this simplified
approach and use a direct trajectory calculator as a contingency strategy when
the planning phase fails to find a solution. Table 3 gives a breakdown of the
number of events, processes and actions in a domain, and the number of objects
in a typical level. As noted above, Angry Birds is an atypical domain as the
agent only has one action (release bird).

Level | Objects | Actions | Events | Processes
---|---|---|---|---
22 | 7-11 | 5(1) | 19-114(17) | 3(3)
25 | 8-12 | 5(1) | 20-88(17) | 3(3)
36 | 8-12 | 5(1) | 30-114(17) | 3(3)
45 | 9-13 | 5(1) | 37-162(17) | 3(3)
46 | 11-15 | 9(1) | 59-199(17) | 5(3)
53 | 10-14 | 5(1) | 28-163(17) | 3(3)
54 | 8-12 | 5(1) | 30-114(17) | 3(3)
57 | 11-15 | 5(1) | 65-220(17) | 3(3)
55 | 29-138 | 13(1) | 430-13772(17) | 7(3)
Table 3: Number of objects and happenings in the grounded domain (and in the
domain template). The main computational load is checking whether events have
occurred.

## Solving Angry Birds with a PDDL+ Planner

In this section, we describe how we used the PDDL+ modeling of Angry Birds
specified above to design an Angry Birds playing agent called _Hydra_. We
detail how Hydra interacts with the game, the PDDL+ planner it uses, and
search techniques we implemented to improve its efficiency.

### From Angry Birds to PDDL+

The game playing API represents an Angry Birds level as a list of labeled
objects and their locations. Hydra supplements this information with
background knowledge of the game and its objects, such as the mass of the
different birds. Then, it automatically translates all the collective relevant
information of the current level into a PDDL+ problem file.

Next, it uses a PDDL+ planner to generate a plan that kills at least one pig
using one bird. If no plan is found in 30s, we execute a default non-planning
action, namely a direct shot at a random pig. If a plan has been found, Hydra
performs the plan. After performing an action, Hydra observes the updated
state of the game, again via the game playing API. Then, it generates a
corresponding PDDL+ problem as before, and generates a new plan from the
current state. This continues until the level ends, either passing the level
or failing to do so before exhausting all the available birds.

#### PDDL+ Planners

Any state-of-the-art PDDL+ planners can be used for solving the Angry Birds
game, such as, SMTPlan+ (Cashmore et al. 2016), UPMurphi (Della Penna et al.
2009), DiNo (Piotrowski et al. 2016) and ENHSP (Scala et al. 2016). However,
during experimentation we faced several challenges due to costly domain-
independent heuristics and lack of ability to incorporate other heuristics for
improved results in the game. Thus, we used our own domain-independent PDDL+
planner in Hydra, which is written in Python and has been used successfully
incorporated for several other domains. Now we will discuss some of the search
algorithms that are supported by our planner and the two heuristics that we
evaluated for the domain.

### Search and Heuristics

Our domain-independent PDDL+ planner supports several search algorithms
including depth-first search (DFS), breadth-first search (BFS), and Greedy
Best-First Search (GBFS)111Our planner also support A* (Hart, Nilsson, and
Raphael 1968), but since we do not aim to find optimal solutions, there is no
benefit in using A*. Nevertheless, we performed an initial evaluation of using
A* in Hydra, but the results were subpar.. GBFS requires heuristic function to
evaluate explored states. We propose two heuristics for our domain.

#### The Score Heuristic

While the agent aims to _solve_ levels in the game, i.e., kill all pigs, the
game itself associates a score to each state. Specifically, destroying pigs
and blocks increases the score of the current level. The score heuristic,
denoted $H_{S}$, uses exactly this score as a heuristic. The agent needs to
analyze the opportune angles at which to release the slingshot and evaluate
distant effects. The heuristic is calculated based on the number of pigs and
blocks the bird will hit multiplied by a constant value.

#### The Proximity Heuristic

The second heuristic we used calculates the active bird’s proximity to the
nearest pig, as a function of flight time and direction. This heuristic,
denoted $H_{P}$, is computed as follows.

$\displaystyle\vec{d}=\vec{x}_{target}-\vec{x}_{bird}$ vector from bird to pig
(2) $\displaystyle\hat{d}=\frac{\vec{d}}{\lVert\vec{d}\rVert}$ unit vector of
same (3) $\displaystyle
H_{P}(s)=\frac{\rVert\vec{d}\rVert}{\vec{v}_{bird}\cdot\hat{d}}$ flight time
to pig (4)

Where $\vec{x}_{target},\vec{x}_{bird}$ and $\vec{v}_{bird}$ are the
coordinate vectors of the pig, the bird and the velocity of the bird
respectively, $\lVert\vec{x}\rVert$ is the Euclidean norm, and $\cdot$ is the
inner product. Using this heuristic results in giving preference to states
where the bird is not only closer to a pig, but also heading towards it,
reducing the effort spent on distant trajectories. As a consequence, this
heuristic also gives preference to low trajectories. The heuristic is
admissible as the straight-line distance to the target underestimating the
flight time.

### Helpful States

Orthogonal to the choice of heuristics, our system incorporates domain-
specific strategies into the search process. This mechanism, that we refer to
as _helpful states_ , is based on the “helpful operators” mechanism from
classical planning (Richter and Helmert 2009). Helpful operators in classical
planning is a mechanism that prioritizes the use of a subset of the operators
while planning. Porting this mechanism to our case requires some adaptation,
since the number of actions in our domain is very limited.

Our _helpful states_ mechanism accepts some definition of what a preferred
state is, and then prioritize expanding such states. In our implementation we
defined a preferred state as one in which the active bird is on a trajectory
that is expected to hit a pig or a TNT block. Checking if the active bird is
on a trajectory to hit an object is done by using the equation for ballistic
motion to create a function $y_{t}=f_{i}(x_{t})$ for each potential target $i$
(pig or TNT block). Given a state, we can calculate whether the bird’s $y$ is
close to the desired value, and if so, mark it as a preferred state. The
search strategy can also be understood as doing Monte-Carlo rollout (Chaslot
et al. 2008) only on a smaller set of states, instead of performing a complete
rollout on all the search states.

Throughout the search we maintain two open lists, one that only includes
preferred states, and one that includes all states. During the search, we
alternate between expanding states from the open list with the preferred
states and the regular open list. To ensure that preferred states are reached
early, we also mark all states where the bird has not yet been launched as
preferred. This ensures the planner reasons with the full range of possible
launch angles. Note that this “preferred states” mechanism is fairly general,
and allows different definitions of preferred states, and multiple open lists
with different priorities. This is a topic for future research.

## Experiments

|  |
---|---|---
22: A level requiring a single low-trajectory shot. |  25: A level requiring a single high-trajectory shot. |  36: In this level, the TNT is too far away to kill the pig.
|  |
45: To kill both pigs with one bird, must hit TNT. |  46: the TNT is a decoy, must hit each pig with a bird. |  53: Must hit TNT to win, but it is a difficult shot.
|  |
54: Bird must bounce off far platform wall to hit pig. |  57: TNT provides a solution to an ‘impenetrable’ box. |  55: A complex level.

Figure 8: A screenshot example of each type of level layout

We conducted a set of experiments to evaluate the efficacy of various
planning-related design decisions explored in this paper to develop Hydra.
Particularly, in experiment 1, we studied how different search perform when
playing Angry Birds using the PDDL+ model introduced in this paper. Next, in
experiment 2, we measured the impact of different heuristics have on
controlling the search. These experiments were conducting on a benchmark set
of Angry Birds problems. We also present results from three other domain-
specific Angry Birds agents on the benchmark set as a way to situate Hydra’s
performance.

#### Benchmark Levels

Our benchmark set of problems contains _simple levels_ that contain a single
bird and a relatively small number of other objects, and _complex levels_ that
contain multiple birds and more than 50 objects. The evaluation is designed to
test the agents in situations requiring varied strategies. The simple levels
require high accuracy shots while the complex levels require a high level of
physics reasoning about interaction with blocks, but can often be passed by
simple curated rules.

In total, we experimented with 9 types of levels, which are numbered 22, 25,
36, 45, 46, 53, 54, 57, and 55. Every evaluated agent attempted to solve 25
randomly generated levels of each type. An example layout for each type of
level can be found in figure 8. Note that the first 8 level types are simple
levels and the last one is a complex level.

#### Competing Agents

We selected a baseline agent by ANU and two former champions of the IJCAI
competition.

  * •

ANU Baseline (Stephenson and Renz 2017) is a baseline agent developed by the
IJCAI AI Birds competition organizers. This baseline agent targets a randomly
selected pig with each active bird, and aims to shoot at the selected pig
while disregarding all other objects in the scene.

  * •

DataLab (Borovicka, Spetlik, and Rymes 2014) is the 2014 & 2015 champion of
the IJCAI Angry Birds competition. DataLab’s selects from 4 predefined
strategies (destroy pigs, target TNT, target round blocks, destroy
structures), based on predicted utilities.

  * •

Eagle’s Wing (Wang 2017) is the 2016-2018 champion of the ICJAI Angry Birds
competition. Eagle’s Wing selects from 5 predefined strategies (pig-shooter,
TNT, most blocks, high round objects, and bottom building blocks), based on
predicted utility. Eagle’s Wing supplements these strategies with a trained
XGBoost (Chen and Guestrin 2016) model to optimize its performance.

  * •

Hydra: We evaluated different variations of Hydra that use different search
strategies, breadth-first search (Hy.,BFS), depth-first search (Hy.,DFS), and
greedy-best first search (Hy.,GBFS). Three variations of Hy.,GBFS were tested
each leveraging different heuristics discussed in this paper: score heuristic
($H_{S}$), proximity heuristic ($H_{P}$), and helpful states ($HS$).

All agents are designed to work with the ANU Science Birds framework (Renz et
al. 2015).

#### Results 1

| Levels
---|---
Agent |  22 |  25 |  36 |  45 |  46 |  53 |  54 |  57 |  55
Hy., BFS |  21 |  21 |  18 |  25 |  18 |  15 |  11 |  25 |  4
Hy., DFS |  21 |  16 |  16 |  21 |  10 |  13 |  11 |  17 |  2
Hy., GBFS($H_{S}$) |  21 |  21 |  20 |  25 |  14 |  15 |  14 |  25 |  1
Hy., GBFS($H_{P}$) |  21 |  21 |  24 |  25 |  25 |  25 |  16 |  25 |  4
Hy., GBFS($HS$) |  21 |  21 |  25 |  25 |  24 |  25 |  16 |  25 |  9
ANU Baseline |  18 |  0 |  0 |  16 |  0 |  2 |  0 |  11 |  17
Datalab |  15 |  8 |  0 |  24 |  3 |  0 |  0 |  22 |  16
Eaglewings |  2 |  6 |  0 |  23 |  2 |  0 |  1 |  21 |  15
Table 4: Evaluation results for each agent, reporting on number of levels
passed.

First, we measure the number of levels _passed_ by each agent. This metric
ensure that we study not only that the agents can plan but are able to find a
plan that can be executed in the environment for successful goal achievement.
Table 4 shows the number of levels passed from each of the level templates in
our benchmark, by each agent.

The results show several clear trends. First, all Hydra agents outperformed
the other agents in all the simple levels. In some cases, e.g., level type 36,
all GBFS variants of HYDRA passed at least 20 out of 25 levels of this type
while none of the non-Hydra agents were able to pass any levels. On the other
hand, the baseline agents outperformed all Hydra agents in the complex levels.

Hydra’s heuristics guide the planner through the physics of the problem, and
are independent of level structure. As a result, Hydra is accurate and
versatile, and it can reason with any level composition. Hydra achieves best
performance in simple levels requiring highly accurate shots.

In contrast, non-Hydra agents were developed for playing human-designed
levels, which are readily solved by a small number of fixed strategies focused
on destroying complex structures. Low-accuracy strategies fare well on
densely-populated levels (i.e., level 55 in the evaluation). Agents developed
for levels with hand-crafted choices of strategy, and requiring only low-
accuracy shots, do not perform well when the underlying assumptions for those
strategies are broken (i.e., randomly generated level compositions).

Level | 22 | 25 | 36 | 45 | 46 | 53 | 54 | 57 | 55
---|---|---|---|---|---|---|---|---|---
Nodes/s | 1,610 | 1,295 | 930 | 624 | 732 | 747 | 1,202 | 634 | 36

Table 5: Rate of node expansion, in nodes per second, by Hydra using GBFS with
helpful states (GBFS (HS)), for each of level type in our benchmark (top row).

Hydra does not perform well in complex levels since the sheer number of
objects present require higher computational resources to reason about them.
However, this suggests that smarter search and heuristics may allow our agent
to also solve the complex levels - an extension we are considering for future
work. Consider Table 5 that lists the number of nodes expanded per second when
using Hydra with GBFS and the helpful states technique. As can be seen, node
generation in the complex level type (type 55) is much slower, going from
1,610 nodes/s in simple level type 22 to only 36 nodes/s for the complex level
type 55. Solving complex levels with domain-independent planning requires
stronger heuristics and perhaps other search techniques.

|  |  | GBFS
---|---|---|---
Level | BFS | DFS | $H_{S}$ | $H_{P}$ | $HS$
22 | 11,468 | 20,000 | 11,495 | 10,081 | 8,445
25 | 36,681 | 46,108 | 35,432 | 21,382 | 9,668
36 | 32,111 | 42,482 | 32,753 | 18,383 | 1,970
45 | 5,257 | 21,341 | 5,279 | 2,264 | 1,826
46 | 31,517 | 40,871 | 31,317 | 15,485 | 2,674
53 | 26,340 | 36,199 | 29,366 | 14,909 | 2,760
54 | 34,564 | 42,100 | 35,115 | 23,578 | 15,526
57 | 6,126 | 17,974 | 7,298 | 2,378 | 1,753
55 | 1,208 | 1,489 | 1,454 | 970 | 707
Table 6: Average total number of nodes expanded by the different Hydra agents
for different level types.

#### Results 2

Next, we study how various heuristics implemented in Hydra impact plan search.
To measure search efficacy, we recorded the number of nodes expanded during
planning. Table 6 shows the average number of nodes expanded by each of the
Hydra agents we considered. Here we see that all of our implemented heuristics
are effective as they reduce the number of nodes expanded during search.

While intuitively appealing, the score heuristic ($H_{S}$) turns out to be
ineffective since the PDDL+ model of the game has only one action (release the
slingshot at an opportune angle), and the action has distant effects. The
score is not a direct consequence of the action, but rather it is computed far
in the future once the bird collides with other objects.

Helpful states mechanism is the most effective; in some cases it expands an
order fewer nodes before finding a solution (see the row corresponding to
level 53 for example). Returning to the results in Table 4, we see that search
efficiency significantly impacts Hydra’s performance. GBFS with various
heuristics can solve more levels that BFS and DFS. Between the different
search enhancements, the helpful states technique is clearly beneficial in
most cases.

## Conclusion

In this work, we presented a novel approach based on domain-independent
planning for Angry Birds, a popular mobile games and a challenging AI testbed
(Renz et al. 2019). To capture the complex dynamics of the system, we modeled
the game using PDDL+, a rich planning language developed for hybrid systems
supporting exogenous activity in planning domains via discrete events and
continuous processes. We discuss how expressiveness of PDDL+ can capture
complex, physics-based dynamics of Angry Birds. We implemented this PDDL+
model within a game playing planning agent, and explored several domain-
specific search enhancements, such as novel heuristics and the “helpful-
states” search strategy. Such improvements were useful to efficiently solve
the challenging Angry Birds domain, despite a PDDL+ model focused on reducing
system complexity and only incorporating a single action. Our evaluation
showed promising results, with our agent able to pass more levels than
multiple champions of the Angry Birds AI competition, in most types of levels
we considered. It highlights that domain-independent planning and
combinatorial search is a viable approach for solving this AI benchmark.

Despite promising results, our current agent can struggle with complex levels
composed of several objects on the scene. Additionally, our implementation
does not consider the full range of the Angry Birds game, ignoring birds with
special abilities and various types of pigs. Future work will address both
limitations: develop better heuristics and search algorithms and extending our
domain model to include birds’ special powers. Future work will also compare
Hydra with emerging AI agents such as BamBirds and Agent X(Lutalo 2022), which
improve shot selection by better balancing between exploration and
exploitation.

Our work demonstrates that PDDL+ is useful for capturing the dynamics of
challenging, real-world domains. However, our experience with Angry Birds
demonstrates that solving planning problems in rich PDDL+ domains is a
challenge for existing planners. The development of solutions for PDDL+ by the
community has been slow. One of the reasons for this is that PDDL+ benchmark
domains have been largely re-purposed from existing domains, which result in
problems that could be addressed by PDDL2.1 techniques. This work also serves
as a call-to-arms for research into PDDL+ planners and heuristics, both
domain-independent and model-specific.

## Acknowledgements

This work was supported by the DARPA SAIL-ON program under contract
HR001120C0040. The views and conclusions in this document are those of the
authors and should not be interpreted as representing the official policies,
either expressly or implied, of the Defense Advanced Research Projects Agency
or the U.S. Government.

## References

  * Alur et al. (1995) Alur, R.; Courcoubetis, C.; Halbwachs, N.; Henzinger, T.; Ho, P.; Nicolin, X.; Olivero, A.; Sifakis, J.; and Yovine, S. 1995.  The Algorithmic Analysis of Hybrid Systems.  _Theoretical Computer Science_ , 138: 3–34.
  * Borovicka, Spetlik, and Rymes (2014) Borovicka, T.; Spetlik, R.; and Rymes, K. 2014.  DataLab Birds Angry Birds AI.  _Czech Technical University in Prague, Prague, Czech Republic, Tech. Rep. of the Datal. ab_.
  * Cashmore et al. (2016) Cashmore, M.; Fox, M.; Long, D.; and Magazzeni, D. 2016.  A Compilation of the Full PDDL+ Language into SMT.  In _ICAPS_ , 583–591.
  * Chaslot et al. (2008) Chaslot, G. M. J.; Winands, M. H.; Herik, H. J. v. d.; Uiterwijk, J. W.; and Bouzy, B. 2008.  Progressive strategies for Monte-Carlo tree search.  _New Mathematics and Natural Computation_ , 4(03): 343–357.
  * Chen and Guestrin (2016) Chen, T.; and Guestrin, C. 2016.  Xgboost: A scalable tree boosting system.  In _ACM International Conference on Knowledge Discovery and Data Mining (SIGKDD)_ , 785–794.
  * Della Penna et al. (2009) Della Penna, G.; Magazzeni, D.; Mercorio, F.; and Intrigila, B. 2009.  UPMurphi: a tool for universal planning on PDDL+ problems.  In _Nineteenth International Conference on Automated Planning and Scheduling_.
  * Edelkamp and Hoffmann (2004) Edelkamp, S.; and Hoffmann, J. 2004.  PDDL2. 2: The language for the classical part of the 4th international planning competition.  _4th International Planning Competition (IPC’04), at ICAPS’04_.
  * Fikes and Nilsson (1971) Fikes, R. E.; and Nilsson, N. J. 1971.  STRIPS: A new approach to the application of theorem proving to problem solving.  _Artificial intelligence_ , 2(3-4): 189–208.
  * Fox and Long (2003) Fox, M.; and Long, D. 2003.  PDDL2.1: An Extension to PDDL for Expressing Temporal Planning Domains.  _Journal of Artificial Intelligence Research_ , 20: 61–124.
  * Fox and Long (2006) Fox, M.; and Long, D. 2006.  Modelling mixed discrete-continuous domains for planning.  _Journal of Artificial Intelligence Research_ , 27: 235–297.
  * Hart, Nilsson, and Raphael (1968) Hart, P. E.; Nilsson, N. J.; and Raphael, B. 1968.  A Formal Basis for the Heuristic Determination of Minimum Cost Paths.  _Systems Science and Cybernetics, IEEE Transactions on_ , 4(2): 100–107.
  * Henzinger (2000) Henzinger, T. A. 2000.  The theory of hybrid automata.  In _Verification of digital and hybrid systems_ , 265–292. Springer.
  * Lutalo (2022) Lutalo, D. 2022.  Agent X: Improving Exploration vs Exploitation in the State of the Art Angry Birds AI.  In _2022 IEEE Conference on Games (CoG)_ , 330–337. IEEE.
  * McDermott et al. (1998) McDermott, D.; Ghallab, M.; Howe, A.; Knoblock, C.; Ram, A.; Veloso, M.; Weld, D.; and Wilkins, D. 1998.  PDDL - The Planning Domain Definition Language.  _Technical Report_.
  * McDermott (2003) McDermott, D. V. 2003.  Reasoning about Autonomous Processes in an Estimated-Regression Planner.  In _International Conference on Automated Planning and Scheduling_.
  * Piotrowski et al. (2016) Piotrowski, W.; Fox, M.; Long, D.; Magazzeni, D.; and Mercorio, F. 2016.  Heuristic Planning for PDDL+ Domains.  In _IJCAI_ , 3213–3219.
  * Ramirez et al. (2017) Ramirez, M.; Scala, E.; Haslum, P.; and Thiebaux, S. 2017.  Numerical integration and dynamic discretization in heuristic search planning over hybrid domains.  _arXiv preprint arXiv:1703.04232_.
  * Renz et al. (2019) Renz, J.; Ge, X.; Stephenson, M.; and Zhang, P. 2019.  Ai meets angry birds.  _Nature Machine Intelligence_ , 1(7): 328–328.
  * Renz et al. (2015) Renz, J.; et al. 2015.  AIBIRDS: The Angry Birds Artificial Intelligence Competition.  In _AAAI_ , 4326–4327.
  * Richter and Helmert (2009) Richter, S.; and Helmert, M. 2009.  Preferred operators and deferred evaluation in satisficing planning.  In _International Conference on Automated Planning and Scheduling (ICAPS)_ , volume 19, 273–280.
  * Scala et al. (2016) Scala, E.; Haslum, P.; Thiébaux, S.; and Ramirez, M. 2016.  Interval-based relaxation for general numeric planning.  In _Proceedings of the Twenty-second European Conference on Artificial Intelligence_ , 655–663.
  * Stephenson and Renz (2017) Stephenson, M.; and Renz, J. 2017.  Creating a Hyper-Agent for Solving Angry Birds Levels.  In _AIIDE_ , 234–240.
  * Stephenson, Renz, and Ge (2020) Stephenson, M.; Renz, J.; and Ge, X. 2020.  The computational complexity of Angry Birds.  _Artificial Intelligence_ , 280: 103232.
  * Wang (2017) Wang, T. J. 2017.  Description for Eagle’s Wing, an AI Angry Bird Playing Agent.  _University of Waterloo_.


# Analytic solutions for irregular diffusion equations with concentration
dependent diffusion coefficients

I.F. Barna<EMAIL_ADDRESS>L. Mátyás  Hungarian Research
Network, Wigner Research Centre for Physics, Konkoly–Thege Miklós út 29-33,
H-1121 Budapest, Hungary  Sapientia University, Department of Bioengineering,
Libertătii sq. 1, 530104 Miercurea Ciuc, Romania

###### Abstract

We investigate diffusion equations which have concentration dependent
diffusion coefficients with physically two relevant Ansätze, the self-similar
and the traveling wave Ansatz. We found that for power-law concentration
dependence some of the results can be expressed with a general analytic
implicit formulas for both trial functions. For the self-similar case some of
the solutions can be given with a formula containing the hypergeometric
function. For the traveling wave case different analytic formulas are given
for different exponents. For some physically reasonable parameter sets the
direct solutions are given and analyzed in details.

††journal: Communication in Nonlinear Science and Numerical Simulation

##  1 Introduction

The simplest mass or energy transfer is diffusion or heat conduction phenomena
which is one of the most important mass transport phenomena in nature
therefore it has an enormous interest in science and engineering as well. The
corresponding literature has grown enormously over the past 200 years. Without
completeness we just mention some recent textbooks [1, 2, 3, 4].

The simplest diffusion process is the regular one, which is formulated with
well-known parabolic partial differential equation (PDE) in the form of

$\frac{\partial C(x,t)}{\partial t}=D\cdot\Delta C(x,t),$ (1)

where $C(x,t)$ is the concentration an D is the diffusion coefficient which is
a positive real constant and $\Delta$ represents the Laplace differential
operator in arbitrary dimensions in arbitrary coordinate system. Certain
boundary conditions belong to equation (1). There are definite solutions for
finite systems, which often are related to engineering applications [5, 6].

There are numerous works done in the last decades to find analytic solutions
beyond the well-known Gaussian and error functions. The best known is the work
of Bluman and Cole [7] from 1969 who gave an analysis based on a general
symmetry analysis method giving numerous analytic solutions, some of them are
expressible with Gaussian or error functions. As we see they did not present
any kind of solutions which looks similar to ours (and which we present in the
following).

In case of infinite horizon the fundamental solution is the Gaussian function
which has application in different areas of science. In the last year we found
additional, physically relevant solutions with the help of the self-similar
Ansatz [14, 15] which are a logical generalization of the fundamental
solution. These solutions are the multiplication of the Gaussian function and
the Kummer’s M and Kumer’s U functions which have an additional free
parameter, due to this parameter even decaying and oscillatory solutions can
be given. The applied self-similar Ansatz will be defined later on in this
study in all details. The diffusion coefficient in certain cases may be
considered constant, however there are cases where it may vary [8]. For two
dimensional models important results have been obtained in ref.[9] and for a
diffusive-reactive case in ref.[10]. Environmental aspects of diffusive very
fine particles is discussed in [11].

Diffusion on surfaces is also a significant topic, with possible irregular
features [12]. The chaotic properties of this latter system have been studied
in [13].

The general form of diffusion equation comes from a conservation law and
reads,

$\frac{\partial C}{\partial t}=\frac{\partial}{\partial
x}\left(D\frac{\partial C}{\partial x}\right).$ (2)

In case the diffusion coefficient $D$ depends on parameter $C$, then in
general it will also depend on $x$. So the case $D(C[x,t])$ is possible. The
diffusion coefficient may depend on certain physical quantities, and it may
vary depending in which phase the system is: gaseous, fluid or solid phase. In
this sense there is a difference between the dependence of mass diffusion
coefficient and heat diffusion coefficient dependence on the parameter $C$.
The $C$ stands for the density or concentration in case of mass diffusion and
$C$ corresponds to temperature in case of heat diffusion.

The investigation of the regular diffusion is just the first step to
understand diffusion process in general, however there are much more complex
and difficult diffusion processes in nature. The diffusion coefficient can
have additional dependencies like, time, space of even a tricky combination of
both. (We examined the case when the diffusion coefficient depend on the
function $x/t^{\frac{1}{2}}$) In the last years we systematically investigated
of such equations and gave new self-similar solutions together with detailed
numerical investigations as well. We applied explicit, semi-explicit and
implicit numerical schemes to solve numerically the corresponding PDEs and
compared the evaluated solutions to the exact mathematical ones [16, 17, 18].
We found that in most cases the leapfrog-hopscotch method is the most
expedient method to solve PDEs. These kind of diffusion equations with time,
space and "time-and-space" dependent diffusion coefficients do have analytic
solutions with the similar structure the fundamental Gaussian solution
multiplied by the Kummer’s M and Kummer’s U functions or with the Whittaker
type functions. Due to the extra free parameter different type of decaying
solutions are always exist with different asymptotics. Some solutions have
additional oscillations as well.

On this way we may go a bit further and ask the question what are the
properties of the diffusion process when the diffusion coefficient directly
depend on the concentration, in the form of $D(C[x,t])$. These cases cover
certain real non-linear diffusion processes described with highly non-linear
PDEs. In the following we will investigate this question in details.

##  2 Theory and Results

It is evident that diffusion is in general a three dimensional process beyond
Cartesian symmetry, however we limit our analysis to a single Cartesian space
dependent equation. In case we consider the equation of heat transfer, then we
have for the heat flux $q=-\kappa T_{x}$, where $\kappa$ is the thermal
conductivity. If $\kappa$ depends on temperature we have for the energy
balance equation,

$\rho c_{p}\frac{\partial T}{\partial t}=\frac{\partial}{\partial
x}\left(\kappa(T)\cdot\frac{\partial T}{\partial x}\right).$ (3)

The function of $\kappa(C[x,t])$ in principle can be any kind of continuous
function with existing first derivative with respect to the concentration
$C(x,t)$.

There are numerous studies available in which various the non-linear diffusion
(or heat conduction) equations were solved with different methods and means
and sometimes analytic solutions are given. Without completeness we cite the
most relevant studies from the last seven decades.

We should start our reference with the work of Fujita [19] from 1952 who gave
analytical solutions when the diffusion coefficient has the form of
$D(C)=D(0)/(1-\lambda C)$.

Later Pattle [20] in 1959 gave solutions with compact support for diffusion
from an instantaneous point source in one, two, or three dimensions. Philip
[21] derived exact solutions for the non-linear diffusion equation when the
concentration has the forms of $D(C)=D_{0}/(1-\lambda C),D(C)=D_{0}/(1-\lambda
C)^{2}$ and $D(C)=D_{0}/(1+2aC+bC^{2})$.

Boyer [22] used a special kind of general self-similar Ansatz with the form of
$T(x,t)=U(t)Y(r/R(t))$ and solved the non-linear heat conduction equation when
the thermal conductivity was given in the form of $\kappa(T)=(\alpha+\beta
T+\gamma T^{2})^{-1}$.

In 1964 Bankoff [23] investigated heat conduction or diffusion with change of
phase and presented numerous methods how the solutions can be derived with
series expansion or integral methods.

Knight and Philip [24] gave solution for the case of $D(C)=a(b-C)^{-2}$ with
the help of the linearisation of the equation.

Tuck [25] solved the diffusion equation for $D(C)=kC^{n}$ with the constant
source boundary condition with the self-similar solution and presented results
with compact supports.

Munier et al. [26] presented that the self-similar and the partially invariant
solutions are identical and introduced the theory of homology with new type of
solutions. These are related through the Bäcklund transformation. They
acknowledge in the paper, that for $n=-1$ is an exceptional case, where they
met singularities. In this paper we discuss the case $n=-1$, and we give an
explicit continuous solution.

King [27] solved a cylindrical symmetric non-linear diffusion equation with
the self-similar Ansatz (which is very similar to our) and found solutions
which can be expressed with the Airy functions.

Sadighi and Ganji [28] applied the variational iteration method and presented
analytic results in form of final polynomials.

Hayek [29] presented an exact solution for a nonlinear diffusion equation in a
radially symmetric n-dimensional case in inhomogeneous medium with the help of
the self-similar Ansatz (Eq.6) which we also apply.

Kosov and Semenov [30] derived new radially symmetric exact solutions of the
multidimensional nonlinear diffusion equation, which can be expressed in terms
of elementary functions, Jacobi elliptic functions, Bessel functions, the
exponential integral and the Lambert-W function.

Luckily to us, it must be strongly emphasized that we cannot identify
solutions in the above cited literature which is exactly the same as our Eq.
(16). Our personal problem with some of the mentioned studies that after doing
the symmetry reduction the mathematical properties of derived ODEs are not
investigated and the parameter dependence of the solutions are seldom
discussed in details. We think that is the crucial point of the whole
investigation and when that is properly done then much larger part of the
scientific community (including applied physicists, chemist or various
engineers) can directly profit from these solutions. We think that clear-cut
well-explained analytic formulas still have high relevance in the
understanding of physics.

As an outlook we mention some analytic studies of some non-linear reaction
diffusion equations which are non-linear diffusion equations with one or more
extra source terms.

An exhaustive group classification of such equation were evaluated by
Dorodnitryn [31] in 1982.

The non-classical symmetry reduction was done by Arrigo et al. [32] in 1994.

Vijayakumar [33] presented a study in which he investigated the generalized
diffusion equations (the Fisher, Newell-Whitehead and Fitzhugh-Nagumo
equations) via the isovector approach and showed analytic results as well.

Cherniha and Serov [34] investigated the more general nonlinear diffusion
equations with convection term with the Lie and non-Lie symmetry methods and
presented additional solutions. However non clear-cut formulas are given with
parameter studies.

Reductions and symmetries properties for a generalized Fisher equation with a
diffusion term dependent on density and space was investigated by Chulian et
al. [35].

Liu [36] gave a generalized symmetry classification, and gave the integrable
properties with exact solutions for some nonlinear reaction-diffusion
equations.

Qu et al. [37] applied the conditional Lie - Bäcklund symmetries with
differential constraints and presented explicit solutions for a class of
nonlinear reaction-diffusion equations.

It is important to mention here that this is just the simplest (the
phenomenological) way to introduce non-linearity into the heat conduction
equation. We still apply the Fourier law, where the heat flux is equal to the
temperature gradient times thermal conductivity $q=-\kappa T_{x}$ which has
now got an extra temperature dependence $\kappa(T)$. There is a mathematically
more rigorous method to derive non-linear heat conduction equations which go
beyond the Fourier law. The heat flux should be approximated with the higher
temporal derivatives of the temperature. If the second term is considered we
arrive to the Cattaneo-Vernotte equation [38, 39, 40], more on such kind of
heat conduction equations can be found in the classical work of Gurtin and
Pipkin [41] or in the work of Joseph and Preziosi [42]. Such heat conduction
equations have finite signal propagation velocity properties. An Euler-
Poisson-Darboux type of non-autonomous time-dependent heat conduction equation
was derived by Barna and Kersner [43] which had solutions with a compact
support. To describe heat pulse experiments "beyond the Cattaneo-Vernotte"
models were applied bye Kovács and Ván [44].

Numerical investigation of Eq. (3) was done by [45] applying the implicit
Euler method. In the following we will investigate this equation with the
reduction mechanism applying two physically relevant trial functions, the
self-similar and the traveling wave Ansätze. It is worth to mention here that
this equation is a bit similar to the porous media equation which has the form
of $U_{t}=\Delta(U^{m})$ where $m>1$ and was heavily investigated in former
times [46, 47, 48]. We have to emphasise that if eq. (2) has an extra source
than we arrive to further scientific fields as are the Fisher equation [50,
51], Turing model [52] or Swift - Hohenberg equation [53] for nonlinear
optics. The diffusion equation has a wide range of applicability in science,
which also includes the theory of pricing [54, 55].

After dividing by $\rho c_{p}$ equation (3) can be also written as a diffusion
equation,

$\frac{\partial T}{\partial t}=\frac{\partial}{\partial
x}\left(D_{h}(T)\cdot\frac{\partial T}{\partial x}\right),$ (4)

where $D_{h}=\kappa/(\rho c_{p})$ is the heat diffusion coefficient.

Because the heat and mass diffusion equation has the same PDE form, we simply
consider the general dependence $D=\zeta(C)$ for (2), which corresponds to
$D_{h}=\zeta(T)$ for (4) and we study the following equation:

$\partial_{t}C=\zeta(C)_{C}\cdot C_{x}^{2}+\zeta(C)\cdot C_{xx},$ (5)

where the subscripts mean the partial derivations, respectively. In the
following we study the consequences of such dependence, when the diffusion may
vary in terms of the parameter $C$.

###  2.1 Self-similar Analysis

Let’s start with the self-similar Ansatz [48, 49, 56] of the form of

$C(x,t)=t^{-\alpha}f\left(\frac{x}{t^{\beta}}\right)=t^{-\alpha}f(\eta),$ (6)

where $f(\eta)$ is the shape function with the reduced variable $\eta$, the
two self-similar exponents $\alpha$ and $\beta$ are responsible for the decay
and spreading of the solutions if both have non-negative values. In the last
decade we generalized this kind of Ansatz to multiple spatial dimension and
applied it to the Rayleigh-Bénard convection problems [57, 58] or to the
heated boundary layer equations [59].

To make and in-depth analysis the functional form of the concentration
dependent diffusion coefficient has to be defined. We start with the most
evident case, the power law dependence:

$\zeta(C[x,t])=a\cdot C(x,t)^{n}\hskip 14.22636pt\textrm{where
n}\varepsilon\mathbb{R}$ (7)

and the constant $a$ has to role to fix the dimension. (For simplicity we fix
its numerical value to unity.)

If we insert this equation into (5) we get the following equation

$\rho c_{p}\left(-\alpha
t^{-\alpha-1}f-\beta\,t^{-\alpha-1}f^{\prime}\right)=nt^{-\alpha
n-\alpha-2\beta}f^{\prime 2}+t^{-\alpha n-\alpha-2\beta}f^{\prime\prime},$ (8)

where prime means derivation in respect to $\eta$.

After simplification with $t^{-\alpha}$ one arrives to:

$\rho c_{p}\left(-\alpha t^{-1}f-\beta\,t^{-1}f^{\prime}\right)=nt^{-\alpha
n-2\beta}f^{\prime 2}+t^{-\alpha n-2\beta}f^{\prime\prime}.$ (9)

Both sides of the equation has the same decay in time if

$-1=-n\alpha-2\beta,$ (10)

or $\alpha=(1-2\beta)/n$.

After simplification with $t^{-1}$, and inserting the value of $\alpha$ in the
equation, we arrive to the ordinary differential equation (ODE) of

$\rho c_{p}\left(-\left[\frac{1-2\beta}{n}\right]f(\eta)-\beta\eta
f(\eta)^{\prime}\right)=nf^{n-1}f^{\prime 2}+f^{n}f^{\prime\prime}.$ (11)

Now we have to make case studies for different $n$s.

###  2.2 Case $n=0$

If $n=0$ we have the differential equation which corresponds to the usual
diffusion equation, where the diffusion coefficient is constant. Although this
is considered a relatively known case, here we present two very recent results
for infinite horizon [60]. Beyond the usual Gaussian solution

$C(x,t)=\frac{1}{\sqrt{t}}e^{-\frac{x^{2}}{4Dt}},$ (12)

there are further relative simple solutions. There is a countable set of even
solution relative to the spatial coordinate, the most simple one is (beyond
Gaussian)

$C(x,t)=\frac{1}{t^{\frac{3}{2}}}e^{-\frac{x^{2}}{4Dt}}\left(1-\frac{1}{2D}\frac{x^{2}}{t}\right).$
(13)

There is also a countable set of odd solutions relative to the spatial
coordinate, a simple one is

$C(x,t)=\frac{x}{t^{\frac{5}{2}}}e^{-\frac{x^{2}}{4Dt}}\left(1-\frac{1}{6D}\frac{x^{2}}{t}\right).$
(14)

###  2.3 Case $n\neq 0$

If $n\neq 0$ the equation (11) is more complex. Unlike the regular diffusion
equation, we have now three parameters $\alpha,\beta$ and $n$ and two of them
remains free, $\alpha=(1-2\beta)/n$. The solution has now the direct form of
$C(x,t)=t^{\frac{1-2\beta}{n}}f(x/t^{\beta})$ The obtained ODE (11) has no
general closed form solution for arbitrary $\beta$ and n. We will see that
only for some special fixed $\beta$ and n combinations give us analytic
results.

In the following case $n=-1$ will be considered with special attention,
because corresponds to a characteristic dependence for dilute systems [8, 10].

####  2.3.1  Case $\beta=0$

We may start with the case of $\beta=0$ which give us an implicit solution of

$\int^{f(\eta)}\pm\frac{na^{2n}(2+n)}{\sqrt{-na^{2n}(2+n)(2a^{2+n}\rho
c_{p}-c_{1})}}da-\eta-c_{2}=0.$ (15)

It turned out after some algebra, that if all four parameters of the integral
$n,c_{1},\rho$ and $c_{p}$ are arbitrary rational numbers, there is a definite
solution which can be expressed with the ${}_{2}F_{1}()$ hypergeometric
function [61]

$\displaystyle\left(n(n+2)f(\eta)^{1+2n}\sqrt{1-\frac{2\rho
c_{p}f(\eta)^{2+n}}{c1}}\right)\times$
${}_{2}F_{1}\left(\frac{1}{2},\frac{1+n}{2+n};1+\frac{1+n}{2+n};\frac{2\rho
c_{p}f(\eta)^{2+n}}{c_{1}}\right)\times$
$\displaystyle\left((1+n)\sqrt{-n(2+n)f(\eta)^{2n}[2\rho
c_{p}f(\eta)^{2+n}-c1]}\right)^{-1}-$ $\displaystyle\eta-c_{2}=0.$ (16)

Considering the more special case of $c_{1}=0$ the integral can be given in
closed form for general $\rho c_{p}$ and $n$ values, so the implicit equation
reads as follows:

$\pm\frac{\sqrt{2}(2+n)f(\eta)^{1+2n}}{\sqrt{-nf(\eta)^{3n+2}(2+n)\rho
c_{p}}}-\eta-c_{2}=0\,.$ (17)

As we can see on this result, that in the denominator of the fraction, the
argument of the square root is positive if for certain constraints. One of the
possibilities is if $n$ is a negative number with small absolute value. We
rise this latter equation to the second power and we get

$f(\eta)^{n}=-\frac{n\rho c_{p}}{2(2+n)}(\eta+c_{2})^{2}.$ (18)

For $n=-1$ both sides of the equation are positive

$\frac{1}{f}=\frac{\rho c_{p}(\eta+2)^{2}}{2}.$ (19)

This gives for the $C(x,t)$, in case $\alpha=(1-2\beta)/n=-1$:

$C(x,t)=tf(\eta)=t\cdot\frac{2}{\rho c_{p}(x+c_{2})^{2}}.$ (20)

One may arrive to this result by applying to equation (2), the standard change
of variables $C(x,t)=A(t)\cdot B(t)$, for $n=-1$. This relation fulfils the
equation (3), however it is divergent for large times.

####  2.3.2 Case $\beta=1$

For the second case let’s take $\beta=1$ and $n=-1$ with the solution of

$f(\eta)=\frac{c_{1}^{2}}{-\rho
c_{p}(1+c_{1}\eta)+c_{1}^{2}c_{2}e^{c_{1}\eta}},$ (21)

It is clear that for the case $\rho
c_{p}(1+c_{1}\eta)=c_{1}c_{1}^{2}e^{c_{1}\eta}$ the shape function becomes
singular at one or two points where the linear equation is touches or
intersects the exponential equation a nice example, when all four parameters
have unit values. We exclude such non-physical solutions from our analysis. It
is also clear that such solutions arise when the $\rho$ and $c_{p}$ are much
larger than the initial conditions $c_{1}$ and $c_{2}$.

Figure (1) shows the solutions of Eq. (21) for the shape functions for two
different parameter sets.

Figure 1: The solution of Eq. (21). The used parameter sets $\rho\cdot
c_{p},c_{1},c_{2}$ for the red and green lines are $(1,1,2)$, and $(1,-2,5)$,
respectively.

We analyze the function (21) by evaluating the derivative of it

$\frac{\partial f}{\partial\eta}=\frac{c_{1}^{2}(\rho
c_{p}c_{1}-c_{1}^{3}c_{2}e^{c_{1}\eta})}{(-\rho
c_{p}(1+c_{1}\eta)+c_{1}^{2}c_{2}e^{c_{1}\eta})^{2}}.$ (22)

As one can see, this function has an extreme value if

$\rho c_{p}c_{1}=c_{1}^{3}c_{2}e^{c_{1}\eta}.$ (23)

This means that this extrema will occur at the value of $\eta$

$\eta_{*}=\frac{1}{c_{1}}\ln\frac{\rho c_{p}}{c_{1}^{2}c_{2}}.$ (24)

One can see, that if

$\displaystyle\rho c_{p}$ $\displaystyle>$ $\displaystyle
c_{1}^{2}c_{2},\,\,\,\eta_{*}>0,$ (25) $\displaystyle\rho c_{p}$
$\displaystyle<$ $\displaystyle c_{1}^{2}c_{2},\,\,\,\eta_{*}<0.$ (26)

Figure (2) however presents the $C(x,t)$ total solution in the form of:

$C(x,t)=\frac{1}{t}\left(\frac{c_{1}^{2}}{-c_{p}\rho(1+c_{1}\rho(x/t)+c_{1}^{2}c_{2}e^{\frac{c_{1}x}{t}}}\right).$
(27)

Figure 2: The solution of Eq. (27), the presented $C(x,t)$ function is for
$\alpha=1,\beta=1,\rho=1,c_{p}=0.1,c_{1}=15,c_{2}=0.7,$ parameter set,
respectively.

As a third mathematical case we found a solutions for the parameter pair of
$\beta=-1/2$ and $n=-2$. Unfortunately, the result is a multi-valued implicit
formula with real and complex parts. We tried to tune the parameters
$c_{1},c_{2}$ and $\rho c_{p}$ but cannot found any solution which could be
interpreted physically e.g. has some reasonable asymptotic for infinite time
and space coordinates.

There are exotic but existing real materials which have temperature dependent
heat conduction coefficients. As a first example we may mention magnetically
aligned single wall carbon nanotube films [63] where the heat conduction
coefficient has linear temperature dependence between 50 and 250 K. Our second
example is bulk semiconductor at large temperature gradient. The authors
approximate the heat flux with a sum of higher spatial derivatives of the
temperature

$q=-\kappa_{0}(T)T_{x}-\kappa_{1}(T)T_{xxx}-\kappa_{2}(T)-\kappa_{3}(T)(T_{x})^{3},$
(28)

where the first coefficient is $\kappa_{0}=1/T$ [64]. In our present model we
cannot take into account the higher terms.

###  2.4 Traveling wave analysis

The second physically relevant trial function which we use is the traveling
wave Ansatz in the form of:

$C(x,t)=g(x+ct)=g(\omega),$ (29)

to avoid further misunderstanding we use a different notation for the shape
function which is $g$ and for the reduced variable which is $\omega$ now.

To make and in-depth analysis the functional form of the concentration
dependent diffusion coefficient has to be defined. Let’s try the most evident
case, the power law dependence first:

$\kappa(C[x,t])=a\cdot C(x,t)^{n}\hskip 14.22636pt\textrm{where n,
b}\>\varepsilon\mathbb{R}\backslash 0,$ (30)

$n$ is a free exponent and $a$ is responsible for the proper physical
dimension of the thermal conductivity. (The numerical value of $a$ is set to
unity again.) After the usual algebraic steps we arrive to the ODE of

$\rho c_{p}cg^{\prime}=a\left(ng^{n-1}g^{\prime
2}+g^{n}g^{\prime\prime}\right).$ (31)

With the help of Maple 12 we can derive a general implicit formula which
contains an integral

$\int_{g(\omega)}\frac{aZ^{n}}{c_{1}a+Z\rho c_{p}c}dZ-\omega-c_{2}=0.$ (32)

Luckily, for $n=-1,0$ and $1$ exist closed form solutions. For $n=0$ we get
back the regular diffusion equation with the exponential front solution, which
is nonphysical. For $n=1$ the solution is the sum of the Lambert W function
[61] with the pure argument of $\omega$ plus a function of $\omega$. All
together the solution is divergent at large x and t arguments. For
completeness we mention that in the work of Kosov and Semenov [30] a
completely different solution is presented where an exponential function has
an argument proportional to $[x^{4}/t+\text{Lambert W}(x^{4}/t)]$. We cannot
transform the two solutions into each other with finite algebraic steps, so
our solution is different to [30]. (One may find more about Lambert W function
in [62]. ) Luckily, for $n=-1$ the solutions become simpler and we get

$g(\omega)=\frac{e^{-\frac{\omega+c_{2}}{c_{1}a}}}{c_{1}\rho
c_{p}ce^{-\frac{\omega+c_{2}}{c_{1}a}}-1}.$ (33)

Using the definition of the traveling wave Ansatz the final form of the
concentration reads:

$C(x,t)=\frac{e^{-\frac{(x+ct)+c_{2}}{c_{1}a}}}{c_{1}\rho
c_{p}ce^{-\frac{(x+ct)+c_{2}}{c_{1}a}}-1}.$ (34)

Figure (3) shows the $C(x,t)$ concentration function Eq. (34) for the set of
parameters $c_{2}=0,\rho\cdot c_{p}=1,c=1,c_{1}=10,a=0.1$.

Figure 3: The concentration function $C(x,t)$ of Eq. (34) for the parameter
set $c_{2}=0,\rho\cdot c_{p}=1,c=1,c_{1}=10,a=0.1$.

As a second class of functions we may consider a temporal and spatially
periodic dependence of $\kappa(T[x,t])=b\sin(T[x,t])$ unfortunately we cannot
derive any solution is a reasonable closed form. Additionally we tried the
Lorenzian form of $\kappa(T[x,t])=\frac{a}{1+T(x,t)^{2}}$ and the exponential
form of $\kappa(T[x,t])=b\cdot Exp(-T[x,t])$ in vain, there are no analytic
closed form available.

##  3 Summary and Outlook

We investigated the highly non-linear diffusion equation where the diffusion
constant (now it is rather a parameter) directly depends on the concentration.
Two type of trial functions were used and different functional form were
analyzed. We found physically relevant analytic solutions which have power law
decays at infinite times. In the future - as a natural generalization - we try
to investigate reactions diffusion equations which are diffusion equations
with extra source terms on the right hand side.

##  4 Acknowledgment

\- One of us (I.F. Barna) was supported by the NKFIH, the Hungarian National
Research Development and Innovation Office.
\- The authors declare no conflict of interest.
\- Both authors contributed equally to every detail of the study.
\- There was no extra external founding.
\- All the evaluated data are available in the manuscript.

## References

  * [1] Ghez G., Diffusion Phenomena; Dover Publication Inc.: Long Island, NY, USA, 2001
  * [2] Vogel G., Adventure Diffusion , Springer, 2019
  * [3] Bird R.B., Stewart W.E. and Lightfoot E.N, Transport Phenomena, 2nd Edition, John Wiley & Sons, Inc., (2002)
  * [4] Lienemann J., Yousefi A., Korvink J.V. Nonlinear Heat Transfer Modeling, In: Benner P. Sorensen D.C. Mehrmann V. (eds) Dimension Reduction of Large-Scale Systems, Lecture Notes in Computational Science and Engineering, 45, (2005) 327 https://link.springer.com/book/10.1007/3-540-27909-1
  * [5] Thambynayagam R.K.M., The Diffusion Handbook: Applied Solutions for Engineers, McGraw-Hill: New York, 2011.
  * [6] Bennett T.D., Transport by Advection and Diffusion  John Wiley & Sons, Hoboken, NJ, 2013
  * [7] Bluman G.W. and Cole J.D. , J. Math. Mech., 18, (1969) 1025
  * [8] Reif F., Fundamentals of Statistical and Thermal Physics, Reissued by Waveland Press, Long Grove, 2009
  * [9] Claus I., Gaspard P., Phys. Rev. E 63, (2001) 036227 https://doi.org/10.1103/PhysRevE.63.036227
  * [10] Mátyás L., Gaspard P., Phys. Rev. E 71, (2005) 036147 https://doi.org/10.1103/PhysRevE.71.036147
  * [11] Salma I., Füri P., Németh Z., Balásházy I., Hofmann W., Farkas Á, Atmospheric Environment 104, (2015) 39 https://doi.org/10.1016/j.atmosenv.2014.12.060
  * [12] Mátyás L., Klages R., Physica D 187, (2004) 165 https://doi.org/10.1016/j.physd.2003.09.008
  * [13] Mátyás L., Barna I.F., Chaos, Solitons & Fractals 44, (2011) 1111 http://dx.doi.org/10.1016/j.chaos.2015.08.002
  * [14] Mátyás L. and Barna I.F. , Rom. J. Phys., 67, (2022) 101 https://rjp.nipne.ro/2022_67_1-2/RomJPhys.67.101.pdf
  * [15] Barna I.F. and Mátyás L. Mathematics, 10, (2022) 3280 https://doi.org/10.3390/math10183281
  * [16] Saleh M., Kovács E., Barna I.F. and Mátyás L. Mathematics 10, (2022) 2813 https://doi.org/10.3390/math10152813
  * [17] Kovács E., Saleh M., Barna I.F. and Mátyás L. Diffusion Fundamentals 35, (2022) 70 https://diffusion.unileipzig.de/pdf/volume35/diff_fund_35(2022)70.pdf
  * [18] Saleh M., Kovács E. and Barna I.F. Algorithms 16, (2023) 183 https://doi.org/10.3390/a16040184
  * [19] Fujita, H. Textile Research Journal, 22, (1952) 823 https://doi.org/10.1177/004051755202201209
  * [20] Pattle. R. E., Quart. Jouro. Mecb. and Applied Math., 12, (1959) 401 https://doi.org/10.1093/qjmam/12.4.407
  * [21] Philip J. R., Australian Journal of Physics, 13, (1960) 1 https://api.semanticscholar.org/CorpusID:117850983
  * [22] Boyer, R. H. Journal of Mathematics and Physics, 40, (1961) 41 https://doi.org/10.1002/sapm196140141
  * [23] Bankoff, S. G., Advances in Chemical Engineering, 5, (1964) 75 https://doi.org/10.1016/s0065-2377(08)60007
  * [24] Knight J. H. and J. R. Philip Journal of Engineering Mathematics, 8, (1974) 219 https://doi.org/10.1007/BF02353364
  * [25] Tuck B, J. Phys. D: Appl. Phys. 9, (1976) 1559 https://doi.org/10.1088/0022-3727/9/11/005
  * [26] Munier, A., Burgan, J. R., Gutierrez, J., Fijalkow, E., and Feix, M. R. SIAM Journal on Applied Mathematics, 40, (1981) 191 https://doi:10.1137/0140017
  * [27] King J.R. J. Phys. A: Math. Gen. 24, (1991) 3213 https://doi.org/10.1088/0305-4470/24/14/010
  * [28] Sadighi A. and Ganji D.D., Computers and Mathematics with Applications, 54, (2007) 1112 https://doi.org/10.1016/j.camwa.2006.12.077
  * [29] Hayek M., Computers and Mathematics with Applications, 68, (2014) 1751 http://doi.org/10.1016/j.camwa.2014.10.015
  * [30] Kosov A. A. and Semenov E. I., Siberian Mathematical Journal, 60, (2019) 93 https://doi.org/10.1134/S0037446619010117
  * [31] Dorodnitsyn V. A. USSR Computational Mathematics and Mathematical Physics 22, (1982) 115. https://doi.org/10.1016/0041-5553(82)90102-1
  * [32] Arrigo D.J., Hill J.M. and Broadbridge P., IMA Journal of Applied Mathematics, 52, (1994) 1 https://doi.org/10.1093/imamat/52.1.1
  * [33] Vijayakumar K., J. Austral. Math. Soc. Set: B 39, (1998) 513 https://doi.org/10.1017/S0334270000007773
  * [34] Cherniha R, and Serov M. Proceedings of the Second International Conference, Symmetry in Nonlinear Mathematical Physics Volume 2, (1997) 444 https://www.slac.stanford.edu/econf/C9707077/papers/symmetry1997.pdf
  * [35] Chulián S., Rosa M. and Gandarias M. L. Journal of Computational and Applied Mathematics 354, (2019) 689 https://doi.org/10.1016/j.cam.2018.11.018
  * [36] Liu H., Commun. Nonlinear Sci. Numer. Simulat. 36, (2016) 21 http://dx.doi.org/10.1016/j.cnsns.2015.11.019
  * [37] Qu G.-Z. Wang G. and Wang M. Results in Physics, 23, (2021) 103971 https://doi.org/10.1016/j.rinp.2021.103971
  * [38] Cattaneo, C. Atti del Seminario matematico e fisico della UniversitÃ di Modenna, 3, (1948) 3
  * [39] Cattaneo, C. Comptes Rendus de l’Académie des Sciences, 247, (1958) 431 https://gallica.bnf.fr/ark:/12148/bpt6k31993/f437.image.langEN
  * [40] Vernotte, P. Comptes Rendus de l’Académie des Sciences, 246, (1965) 3154
  * [41] Gurtin M.E. and Pipkin A.C., Archive for Rational Mechanics and Analysis volume 31, (1968) 113 https://link.springer.com/article/10.1007/BF00281373
  * [42] Joseph, P. P. and Preziosi, L. Reviews of Modern Physics, 60, (1989) 41 https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.61.41
  * [43] Barna I.F. and Kersner R., J. Phys. A: Math. Theor. 43, (2010) 375210 http://iopscience.iop.org/1751-8121/43/37/375210
  * [44] Kovács R. and Ván P. International Journal of Heat and Mass Transfer, 83, (2015) 613 https://doi.org/10.1016/j.ijheatmasstransfer.2014.12.045
  * [45] Filipov S.M. and Faragó I. https://arxiv.org/abs/1811.06337
  * [46] Muskat, M. The Flow of Homogeneous Fluids Through Porous Media. New York: McGraw-Hill. 1937
  * [47] Barenblatt G.I. Prikladna ja Matematika i Mekhanika, 16, (1952) 67
  * [48] Zeldovich, Y.B. and Raizer, Y.P. Physics of Shock Waves and High Temperature Hydrodynamic Phenomena (1st ed.). Academic Press 1966
  * [49] Ames W.F., Nonlinear Partial Differential Equations in Engineering, vol. I,II, Academic Press, New York, 1965.
  * [50] Fisher R.A., Annals of Eugenics 7, (1937) 353 https://doi.org/10.1111/j.1469-1809.1937.tb02153.x
  * [51] Al-Khaled K., J. Comput. Appl. Math. 137, (2001) 245 https://core.ac.uk/download/pdf/82237548.pdf
  * [52] Turing A., Philosophical Transactions of the Royal Society of London B. 237, (1952) 37 https://royalsocietypublishing.org/doi/pdf/10.1098/rstb.1952.0012
  * [53] Arecchi F.T., Boccaletti S., Ramazza PL., Physics Reports 318, (1999) 1 https://doi.org/10.1016/S0370-1573(99)00007-1
  * [54] Mazzoni T., A First Course in Quantitative Finance, Cambridge University Press, Cambridge, 2018
  * [55] Lázár E., Acta Universitatis Sapientiae, Economics and Business 2, (2014) 75 https://intapi.sciendo.com/pdf/10.2478/auseb-2014-0011
  * [56] Sedov L. Similarity and Dimensional Methods in Mechanics, CRC Press 1993.
  * [57] Barna I.F. and Mátyás L. Chaos, Solitons and Fractals 78, (2015) 249 http://dx.doi.org/10.1016/j.chaos.2015.08.002
  * [58] Barna I.F., Pocsai M.A. Lökös S. and Mátyás L. Chaos, Solitons and Fractals 103, (2017) 336 http://dx.doi.org/10.1016/j.chaos.2017.06.024
  * [59] Barna I.F., Bognár G., Mátyás L. and Hriczó K., Journal of Thermal Analysis and Calorimetry 147, (2022) 13625 https://doi.org/10.1007/s10973-022-11574-3
  * [60] Mátyás L. and Barna I.F., Universe 9, (2023) 264 https://doi.org/10.3390/universe9060264
  * [61] Olver F.W.J., Lozier D.W., Boisvert R.F. and Clark C.W. NIST Handbook of Mathematical Functions, Cambridge University Press, 2010 https://dlmf.nist.gov/
  * [62] Mező I., The Lambert W Function Its Generalizations and Applications, CRC Press, 2022
  * [63] Hone J. et all. Appl. Phys. Lett. 77, (2000) 666 https://doi.org/10.1063/1.127079
  * [64] Ezzahri Y., Ordonez-Miranda J. and Joulain K. Int. J. Heat Mass Transf 108, (2017) 1357 http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.01.024