# Weather Sensitive High Spatio-Temporal Resolution Transportation Electric
Load Profiles For Multiple Decarbonization Pathways

Samrat Acharya1, Malini Ghosal1, Travis Thurber1, Casey D. Burleyson1, Yang
Ou2, Allison Campbell1, Gokul Iyer2, Nathalie Voisin13, and Jason Fuller1

{samrat.acharya<EMAIL_ADDRESS>1 Pacific Northwest National
Laboratory, Richland, WA 99354, USA 2 Joint Global Change Research Institute,
Pacific Northwest National Laboratory, College Park, MD 20740, USA 3
University of Washington, Seattle, WA 98195, USA

###### Abstract

Electrification of transport compounded with climate change will transform
hourly load profiles and their response to weather. Power system operators and
EV charging stakeholders require such high-resolution load profiles for their
planning studies. However, such profiles accounting whole transportation
sector is lacking. Thus, we present a novel approach to generating hourly
electric load profiles that considers charging strategies and evolving
sensitivity to temperature. The approach consists of downscaling annual state-
scale sectoral load projections from the multi-sectoral Global Change Analysis
Model (GCAM) into hourly electric load profiles leveraging high resolution
climate and population datasets. Profiles are developed and evaluated at the
Balancing Authority scale, with a 5-year increment until 2050 over the Western
U.S. Interconnect for multiple decarbonization pathways and climate scenarios.
The datasets are readily available for production cost model analysis. Our
open source approach is transferable to other regions.

##  I Introduction

The transportation sector is a major greenhouse gas (GHG) emitter, accounting
for 29% of United States (U.S.) [1] and 27% of global emissions in 2019 [2].
It is the largest GHG emitter in the U.S. and second largest globally. While
the focus has primarily been on electrifying light-duty vehicles (LDVs), the
electrification of medium- and heavy-duty vehicles (MHDVs) and non-road
vehicles (such as trains, aviation, and ships) is gaining momentum [3].
Currently, LDVs have seen significant electrification, with global electric
car sales reaching 16.5 million in 2021, representing 9% of the global car
market [4]. Electric bus and truck penetration is lower at 4% and 0.1%
respectively [4]. However, achieving comprehensive decarbonization requires
electrifying the entire transportation sector, including MHDVs and non-road
vehicles.

Rapid large-scale transportation electrification poses challenges to economics
[5], operations [6], and cybersecurity [7] of power grids and EV chargers. To
address these challenges, accurate and highly granular projections of spatio-
temporal transportation charging load profiles are essential for effective
planning and decision-making. Recent studies [8, 9, 10] have employed data-
driven approaches to project EV charging load profiles. For example, Gaete et
al. [8] developed LDV charging profiles in Germany using data on mobility,
battery characteristics, and charging strategies.

This paper is accepted for publication in 2024 IEEE ISGT NA. The complete
copyright version will be available on IEEE Xplore when the conference
proceedings are published.
---

Borlaug et al. [9] created average daily charging profiles for depot-based
heavy-duty trucks based on truck mobility and charging data. Wang et al. [10]
generated daily charging profiles for MHDVs in California by considering
mobility patterns and future charger capacities. In contrast to these studies
[8, 9, 10], our approach encompasses the entire transportation sector,
providing comprehensive time-series data for charging load profiles.

Figure 1: Energy use of LDVs, MDVs, HDVs, and non-road vehicles in the western
U.S. interconnection for the Business As Usual (BAU) and Net Zero (NZ)
decarbonization pathways in the Global Change Analysis Model.

In this paper we analyze transportation electrification in the western U.S.
interconnection region. Our analysis encompasses two decarbonization pathways:
i) Business As Usual (BAU) and ii) Net-Zero (NZ), modeled using the U.S.
version of the Global Change Analysis Model (GCAM-USA) as depicted in Fig. 1.
GCAM is a widely utilized, economically-driven, multi-sector dynamics model
that simulates the interactions between human and natural systems, playing a
crucial role in assessing global changes and their impacts. Notably, GCAM has
been employed by organizations like the Intergovernmental Panel on Climate
Change (IPCC) [11] and the U.S. Department of Energy. The BAU scenario in GCAM
reflects existing GHG emission reduction policies. The NZ scenario embodies
ambitious plans to achieve a clean U.S. grid by 2035 and net-zero GHG
emissions by 2050. Leveraging these decarbonization pathways, we generate
hourly transportation electric load profiles across Balancing Authorities
(BAs) in the western U.S. interconnection. BAs oversee resource planning,
electricity supply-demand balance, and real-time frequency regulation in
specific geographic areas. For spatial details on BAs in the western U.S.
interconnection, we refer to [12]. To this end, we make the following unique
contributions:

  1. 1.

We develop a spatially-distributed statistical downscaling approach,
leveraging multi-sector dynamics in GCAM, to project annual transportation
energy onto hourly BA-level time-series profiles in the western U.S.
interconnection.

  2. 2.

We conduct a detailed sensitivity analysis of transportation charging profiles
considering climate change scenarios, decarbonization pathways, and charging
strategies.

  3. 3.

We provide open source code at https://doi.org/10.5281/zenodo.7888569[13] to
accelerate the dissemination of the approach and support the community of
practice.

##  II Data and Tools

This section outlines the key data and tools employed to project spatio-
temporal transportation charging load profiles.

###  II-A Global Change Analysis Model (GCAM)

The Global Change Analysis Model (GCAM) [14] is an open-source integrated
assessment model that encompasses energy, economy, agriculture and land use,
water, and climate systems at global spatial scales. For this study GCAM spans
from 2015 to 2100, with 5-year time-steps, and the final calibration year is
2015. We utilize GCAM-USA v6 which focuses on 50 states plus the District of
Columbia and simulates sub-national economy and energy systems. GCAM-USA v6
incorporates socioeconomic drivers, energy transformation trends, and final
energy services at the state level. It utilizes the Shared Socioeconomic
Pathway-2 (SSP2) growth assumptions [15]. GCAM-USA v6 incorporates updated
electricity technology cost assumptions from the National Renewable Energy
Laboratory’s (NREL) 2022 Annual Technology Baseline (ATB) [16]. Transportation
cost and energy intensity assumptions, including those for EVs, are primarily
based on NREL’s Electrification Futures Study [17].

###  II-B Thermodynamic Global Warming Simulations (TGW)

The Thermodynamic Global Warming (TGW) simulations [18, 19] are based on the
Weather Research and Forecasting (WRF) model [20, 21]. WRF is a widely used
numerical weather prediction model that calculates atmospheric variables at
horizontal and vertical grid cells above the Earth’s surface [20, 21]. The TGW
simulations provide hourly and three-hourly meteorological data at a 12 km2
resolution, covering the conterminous U.S. and parts of Canada and Mexico. For
2020-2099, the TGW simulations replay historical weather events under
different levels of global warming based on multiple Representative
Concentration Pathways (RCPs) and Global Climate Models. Warming levels are
derived from average temperature and humidity changes from models that are
“cooler” and “hotter” compared to the multi-model mean. We use four future
climate scenarios obtained by combining the cooler and hotter scenarios with
RCP 4.5 and RCP 8.5.

###  II-C Total ELectricity Loads Model (TELL)

The Total Electricity Loads (TELL) model [22] downscales annual state-level
electricity demand projections from GCAM-USA to a hourly resolution. TELL
utilizes multilayer perceptron models trained on historical hourly electricity
demands and meteorological variations for each BA. The hourly BA-level loads
from TELL are scaled to match the annual state-level loads from GCAM-USA,
providing hourly non-transportation electricity demands for each BA. Combining
these with transportation charging demands allows for a comprehensive analysis
of change in total loads over time and across decarbonization scenarios.

###  II-D EVI-Pro Lite

EVI-Pro Lite [23] is a data-driven tool that projects the aggregated charging
demand of EVs. It uses detailed data on travel patterns, EV attributes, and
charging infrastructure characteristics to generate charging electric load
profiles. The tool is based on advanced PEV simulations trained on real-world
driving data from large U.S. travel databases. Inputs to the tool include
fleet size, discrete daily mean temperature, EV distribution, battery size,
charging access, charging preferences (home vs work), and charging levels. We
acknowledge that future technological advancements (e.g., wireless charging),
changes in charging behavior, and broader EV adoption may challenge EVI-Pro
assumptions based on historical data. Thus, we select parameters to represent
realistic future charging scenarios. We detail our input assumptions and their
relation with GCAM-USA, TGW, and TELL in Section III-A.

###  II-E Fleet DNA Data

The Fleet DNA dataset [24] contains anonymous data on the mobility of
commercial MHDVs in the U.S., including delivery vans, transit buses, and
refuse trucks. We utilize specific data headings, such as ”Vid” (vehicle
identifier), ”start_ts” (vehicle record start time), ”end_ts” (vehicle record
end time), and ”distance_total” (distance traveled by vehicle). To focus on
return-to-base schedules, we exclude mobility instances that are less likely
to be return-to-base. Also, we filter the data to include only vehicles that
return to the depot by midnight and charge during their dwelling time in
depots.

###  II-F Non-Road Vehicles Data

Non-road vehicles data includes enplanements at U.S. commercial airports
published by the U.S. Federal Aviation Administration [25], route miles of
rails, and shipping docks in U.S. published by the U.S. Department of
Transportation [26].

##  III Methodology

Fig. 2 illustrates our methodology for developing spatio-temporal charging
load profiles across BAs in the western U.S. interconnection. Below, we detail
our approach.

Figure 2: Workflow for generating transportation charging load profiles across
BAs in the western U.S. interconnection. Tempr.: Temperature, WRF: Weather
Research and Forecasting Model, TELL: Total ELectricity Loads Model.

###  III-A Developing LDV Charging Load Profiles

Fig. 2 illustrates the generation of LDV charging load profiles. This section
outlines the methodology for translating the state-level annual energy
consumption of electrified LDVs in GCAM-USA into hourly load shapes at the BA
level.

#### Spatial Downscaling

We use county-level resolution data to spatially downscale LDV charging load
profiles. This allows us to incorporate temperature data from the TGW
simulations (Section II-B) and consider the diverse climatic regions within
the BAs. The state-to-county distribution of energy use is determined by
analyzing the distribution of registered EVs across counties using vehicle
distribution data from 2018-2021. Each county’s proportion of the state’s
total EVs is calculated, following the methodology outlined in [5].

#### Temporal Downscaling

The county-level aggregate charging load profiles require various inputs from
different modules. One crucial input for the EVI-Pro module is the fleet size.
While direct fleet size data is not available from GCAM-USA, the annual energy
consumption (in PJ/yr), energy usage per vehicle travel distance (in MJ/kvm),
and average travel distance are utilized to estimate the fleet size. Daily
mean temperature from TGW simulations (Section II-B) is used in EVI-Pro, with
mapping to the nearest permissible discrete values. Additional inputs for EVI-
Pro are outlined in Table I, independent of other modules in this study.

While EVI-Pro [23] does not cover all future charging scenarios, we outline
our key parameter choices in the model.

  * •

Battery EVs have a 250-mile range- the max. battery size in EVI-Pro. Plug-in
hybrid EVs are excluded in this study.

  * •

Predominant use of level-2 chargers (208-240 V) over level-1 chargers (120 V
AC).

  * •

Anticipated decrease in home charging preference from 80% to 60%, influenced
by increasing chargers in public and work. Also, as EV ownership reaches
75-90%, chargers in multi-family residential dwellers increases, reducing home
charging preferences. In 2022, $\approx 63\%$ of houses are single-family
units [27], which also bounds the preference for home charging.

  * •

Study of two charging strategies: _min_delay_ for immediate maximum-speed
charging after arrival and _load_level_ for slow charging during the dwelling
time. By 2035, 30% of LDVs are expected to use load leveling, which is
projected to increase to 70% by 2050. Managed charging with price incentives
and larger battery capacities reduce range anxiety, increasing load-leveling
adoption.

###  III-B Developing MHDV Charging Load Profiles

Fig. 2 illustrates the two-step process for generating MHDV charging load
profiles.

In Step 1, using the MHDV load shape generator in Fig.2, we generate separate
normalized MHDV charging load profiles for MDVs and HDVs. Inputs to the
generator include mobility data for MHDVs (delivery vans, delivery trucks,
school buses, transit buses, bucket trucks, tractors, and refuse trucks)
obtained from Fleet DNA (SectionII-E), MHDV charging strategies, EVCS
capacities ([]kW), unit energy use ([]kWh/miles), fleet size, and number of
sample fleets [9]. Three charging strategies—immediate, delay, and constant
minimum power—are utilized. Immediate charging begins upon depot arrival and
continues until the battery is fully recharged or the next trip commences.
Delayed charging delays charging to ensure full recharge just before the
subsequent trip. Minimum power charging involves charging at a constant
minimum power level throughout the depot dwelling period to guarantee a full
recharge for the next trip. We incorporate the uncertainty of future EV
charging, considering a mix of charging strategies and multiple charging
capacities, by assigning weights to both the charging strategies and charger
capacities. The unit energy use ([]kWh/mi) of vehicles depends on their
weight, which is determined using a survey [28]. To minimize bias in MHDV
profiles, we simulate the MHDV load shape generator (Fig. 2) for a user-
defined fleet size and number of fleet samples, following the approach in [9].
The load shape generator in Fig.2 aggregates Fleet DNA mobility data by month
and generates normalized average daily MHDV charging shapes with a 1-hour
resolution for each month. These daily shapes are extended to create
normalized yearly MHDV charging load profiles. Due to insufficient mobility
data, we assume that the daily MHDV charging shape remains consistent
throughout a month.

In Step 2, we scale the normalized MHDV charging load shapes by the annual
MHDV energy use (in []EJ) projected by GCAM-USA in the western U.S.
interconnection for a given year and decarbonization pathway. The scaled MHDV
charging profiles in the western U.S. interconnection are further downscaled
to BAs by considering the relative penetration of MHDVs across BAs. The
penetration of MDVs in some BAs is reported in [5]. For the remaining BAs we
extrapolate the penetration of MDVs using the penetration of LDVs as a
reference. Additionally, we assume that the relative penetration of MDVs and
HDVs is the same across BAs.

TABLE I: Choice of input parameter for EVI-Pro

PEV | BEV250 | Class | Equal
---|---|---|---
Preference | Home60 | Home access | HA75, HA100
Home power | MostL2 | Home ch. strategy | min_delay, load_leveling
Work power | MostL2 | Work ch. strategy | min_delay

###  III-C Developing Non-Road Vehicles Charging Load Profiles

As Fig.2 shows, we generate constant charging profiles for non-road vehicles
(aviation, rails, and ships). To downscale the state-wise yearly energy use by
non-road vehicles in GCAM-USA to hourly charging power across BAs in the
western U.S. interconnection, we use the ratio of i) airport enplanements for
aviation, ii) route miles travelled in railroads and transits for rails and
trains, and iii) shipping docks for ships, in a given BA to that of the entire
western U.S. interconnection. The energy at the western U.S. interconnection
level is determined by aggregating the state-level energies in GCAM-USA for
the 11 states within the region. In this study, we assume that electrified
rails, ship, and aviation consume constant power over each hour of the year.
We make this assumption due to the lack of data on their charging behavior as
they are not yet widely deployed. Therefore, this study may not reflect
accurate temporal downscaling of non-road vehicles. However, their spatial
downscaling is based on county-level data.

##  IV Case Study and Results

This section examines transportation load profiles in the western U.S.
interconnection for 2035 and 2050, considering different decarbonization
pathways and climate scenarios. It compares the sensitivity of these profiles
and projected loads to the total system load. We present profiles for NZ
decarbonization, highlighting its aggressiveness compared to BAU pathway. Our
open-source code allows for profile generation in 5-year increments until 2050
for both BAU and NZ pathways.

###  IV-A Sensitivity of LDV Load Profiles

Fig. 3 demonstrates the impact of LDV charging strategies, temperature, and
weekdays vs. weekends on the LDV charging profiles. In Fig. 3(a), load-
leveling charging (C2) shows a flat, grid-friendly charging profile compared
to immediate charging (C1). With a 70% adoption rate of C2 by 2050, the peak
load can be reduced by 40% compared to the C1-only charging. Fig. 3(b) shows
the impact of temperature on LDV charging. A temperature of $20^{\circ}$C
requires the least charging, while $-10$, $0$, and $40^{\circ}$C days consume
$28\%$, $21\%$, and $36\%$ more power, respectively. Temperature affects
charging efficiency and energy consumption for the same travel distance. Fig.
3(c) shows LDV charging peak load during weekends is approximately 10% lower
than on weekdays.

(a) (b) (c)

Figure 3: Sensitivity of LDV average daily charging load profiles (00:00 to
24:00 hours in PDT) to (a) charging strategies, (b) temperature, and (c)
weekends and weekdays. In (a), C1 and C2 are immediate (_min_delay_) and load
levelling (_load_level_) charging strategies. The results in (a)-(c) are
simulated for King County, Washington in 2035 using the NZ decarbonization
pathway.

###  IV-B Sensitivity of MHDV Load Profiles

Fig. 4 shows the sensitivity of average daily MHDV charging load profiles to
charging strategies, charger capacities, and number of sample fleets for
stochastic simulation. As Figs. 4(a) and 4(d) show, the choice of charging
strategy significantly impacts charging power’s magnitude and timing,
resulting in varying charging peaks. For instance, as Fig. 4(d) shows, the HDV
peak is 0.67 GW for minimum power charging at 00:00 while it is 2.84 GW for
delay charging at 9:00; representing a 323% increase. To achieve a more
balanced and grid-friendly charging profile, we mix 40% immediate, 10% delay,
and 50% constant minimum power charging strategies as an example. This mix
means majority of MHDVs charge with constant power throughout their dwelling
time or start charging upon arrival at depots and a small portion wait for a
suitable time to charge, such as during periods of cheaper electricity prices,
before their departure. Charger capacities also influence the charging
profiles, with slower charging leading to smaller peaks. For example, as
Fig.4(e) shows, the 50 kW charger exhibits a peak of 0.76 GW at 20:00, while
fast chargers have a peak of 1.08 GW at 19:00, representing a 42% increase. We
anticipate that larger battery capacities will drive the adoption of high
power chargers. Thus, we consider a mix of chargers: 5% 50 kW, 5% 125 kW, 10%
250 kW, 40% 350 kW, and 40% 500 kW. Furthermore, as Figs. 4(c) and 4(f) show,
the charging profiles remain stable across different fleet samples, despite
variations in the number of sample fleets.

(a) (b) (c) (d)

(e) (f)

Figure 4: Sensitivity of MHDV average daily charging profiles (00:00 to 24:00
hours in PDT), (a)-(c): MDV and (d)-(f): HDV, with (a) and (d): charging
strategies, (b) and (e): charger capacities, and (c) and (f): number of sample
fleets for stochastic simulation. In (a) and (d), Imm: Immediate charging,
Del: Delayed charging, Min: Minimum power charging, Ours: 40% Imm + 10% Del +
50% Min strategies. In (b) and (e), Ours: 5% 50 kW + 5% 125 kW + 10% 250 kW +
40% 350 kW + 40% 500 kW chargers. The profiles in (a)-(f) are simulated for
the CISO BA in 2035 using the BAU pathway.

###  IV-C Transportation Load Profiles

Figs. 5(a) and (b) depict the charging load profiles (00:00 to 23:00 UTC hours
discontinuous for each month) for LDVs, MHDVs, rails, aviation, ship, and
total transportation in the western U.S. interconnection in 2035 and 2050
under the NZ decarbonization pathway. Several key dynamics are observed.
Firstly, the charging load profiles differ by vehicle type. Secondly, the
transportation peak load significantly increases from 27.77 GW in 2035 to
59.85 GW in 2050, marking a 115% surge. Additionally, the variation in
charging load experiences a 6% increase, with values of 19.19 GW in 2035 and
20.35 GW in 2050. Thirdly, LDVs contribute prominently to the total
transportation charging load profile. Fourthly, charging load for non-road
transportation modes (aviation, rail, and ship) exhibits an upward trend from
2035 to 2050. Finally, the MHDV charging load slightly decreases in 2050
compared to 2035 (e.g., HDV peak in 2035 is 4.9 GW but 4.1 GW in 2050) due to
the NZ decarbonization pathway’s emphasis on transitioning from MHDVs,
particularly MDVs, to mass transportation options such as rail and ships. This
decline is further driven by the adoption of alternative clean fuel
technologies, such as hydrogen, in later years.

###  IV-D Transportation Electric Load Relative to System Load

Three metrics (M1, M2, and M3) in Fig. 6 assess the spatial and temporal
impact of transportation load on the electric power system. M1 represents the
yearly ratio of total electric energy to transportation electric energy, M2
measures the ratio of transportation load to system load at the system peak,
and M3 measures the ratio of transportation load to system load at the
transportation peak. M1 rapidly increases in 2050 relative to 2035, indicating
steep increase in electrification across all BAs. The results also reveal that
all metrics vary significantly across BAs. For instance, southern California
BAs such as LDWPD and IID have high M1 (e.g., in IID the M1 is 0.50 in 2035
and 0.68 in 2050). This is attributed to the aggressive EV adoption trend in
the region and the predominantly residential system load with a moderate
climate. Similarly, CHPD in Chelan County, Washington, has high M1 due to its
low non-transportation electric load and ambitious EV policy. Notably, the
transportation peak and system peak do not coincide, reducing system stress
from transportation charging.

Figure 5: Average daily transportation charging load profiles (00:00 to 23:00
UTC hours) in the western U.S. interconnection across months in (a) 2035 and
(b) 2050 for the NZ decarbonization pathway and RCP 4.5 climate scenario.

###  IV-E Load Profiles with BAU Decarbonization

In comparison to the BAU pathway (not shown in this paper), the NZ pathway
(Figs. 5 and 6) demonstrates significant differences in charging peak load and
variation. In 2035, NZ has a 4.7% higher charging peak load than BAU (27.77 GW
vs. 26.51 GW), increasing to approximately 65.4% in 2050 (59.85 GW vs. 36.17
GW). NZ also exhibits higher load variation, with about 21.5% more variation
in 2035 (19.19 GW vs. 15.79 GW) and 10.5% more variation in 2050 (20.35 GW vs.
18.41 GW) compared to BAU. Additionally, the transportation-to-system energy
ratio is notably higher in the NZ pathway, particularly in 2050. For example,
in IID, NZ shows a 9.8% higher ratio than BAU in 2050 (0.67 in Fig. 6(b) vs.
0.61).

##  V Conclusion

This paper presents a novel approach for generating spatially-distributed
hourly time-series of transportation charging load profiles. Applied to the
western U.S. interconnection in 2035 and 2050, our analysis reveals that
transportation charging loads can contribute significantly to the system
electric peak, ranging from 2.4% to 56.6% in different BAs, despite accounting
for less than 20% and 17% of total electric load for the 2050 NZ case and the
2050 BAU case, respectively. This variation is influenced by factors such as
load nature, climate zones, and EV adoption trends. Understanding this non-
uniform spatial and temporal impact is crucial for effective decarbonization
and transportation electrification policies. To facilitate further research
and policy analysis, we provide a publicly available dataset and code at
https://doi.org/10.5281/zenodo.7888569. Researchers and decision-makers can
use this resource to develop charging load profiles for other regions,
enabling informed decision-making in the pursuit of sustainable transportation
electrification.

## Acknowledgment

This research was supported by the GODEEEP Investment at Pacific Northwest
National Laboratory (PNNL). PNNL is a multi-program national laboratory
operated for the U.S. Department of Energy (DOE) by Battelle Memorial
Institute under Contract No. DE-AC05-76RL01830. We thank Kate Forrest and
Brian Tarroja at University of California, Irvine, USA for the discussions on
the transportation profiles.

Figure 6: Comparison of electric transportation (trans.) to system load for
the western U.S. interconnection (WI) BAs in (a) 2035 (b) 2050 for the NZ
decarbonization pathway and RCP 4.5 climate scenario.

## References

  * [1] Sources of greenhouse gas emissions. Accessed: Oct 25, 2022. [Online]. Available: https://tinyurl.com/ywjxmx9y
  * [2] Greenhouse gas emissions from energy: Overview. Accessed: Dec 16, 2022. [Online]. Available: https://tinyurl.com/52hk5s32
  * [3] T. T. Mai _et al._ , “Electrification futures study: Scenarios of electric technology adoption and power consumption for the united states,” NREL, Golden, CO (United States), Tech. Rep., 2018.
  * [4] Global EV outlook 2022. Accessed: Jan 05, 2023. [Online]. Available: https://tinyurl.com/4syvj5uc
  * [5] Kintner-Meyer _et al._ , “Electric vehicles at scale–phase I analysis: high EV adoption impacts on the Western US power grid,” Tech. Rep., 2020.
  * [6] J. A. P. Lopes _et al._ , “Integration of electric vehicles in the electric power system,” _Proc. of the IEEE_ , vol. 99, no. 1, pp. 168–183, 2010.
  * [7] S. Acharya _et al._ , “Public plug-in electric vehicles+ grid data: Is a new cyberattack vector viable?” _IEEE Trans Smart Grid_ , vol. 11, no. 6, 2020.
  * [8] C. Gaete-Morales _et al._ , “An open tool for creating battery-electric vehicle time series from empirical data, emobpy,” _Scientific data_ , vol. 8, no. 1, pp. 1–18, 2021.
  * [9] B. Borlaug _et al._ , “Heavy-duty truck electrification and the impacts of depot charging on electricity distribution systems,” _Nature Energy_ , vol. 6, no. 6, pp. 673–682, 2021.
  * [10] B. Wang _et al._ , “Medium and heavy-duty electric vehicle infrastructure \- load operations and deployment (hevi-load) v1.0,” 7 2021. [Online]. Available: https://www.osti.gov/biblio/1836645
  * [11] Climate Change 2022 Mitigation of Climate Change. Accessed: May 8, 2023. [Online]. Available: https://tinyurl.com/cv7pfj6u
  * [12] Western interconnection balancing authorities. Accessed: May 10, 2023. [Online]. Available: https://tinyurl.com/h2wyd4a4
  * [13] S. Acharya, T. B. Thurber, and M. Ghosal, “Transportation Electrification Load Profiles by Balancing Authority in the Western United States for GODEEEP,” May 2023. [Online]. Available: https://doi.org/10.5281/zenodo.7888569
  * [14] GCAM v6 Documentation. Accessed: May 8, 2023. [Online]. Available: https://jgcri.github.io/gcam-doc/toc.html
  * [15] B. C. O’Neill _et al._ , “The roads ahead: Narratives for shared socioeconomic pathways describing world futures in the 21st century,” _Global environmental change_ , vol. 42, pp. 169–180, 2017.
  * [16] Annual energy outlook 2022. Accessed: Feb 19, 2022. [Online]. Available: https://www.eia.gov/outlooks/aeo/
  * [17] P. Jadun _et al._ , “Electrification futures study: End-use electric technology cost and performance projections through 2050,” 2017.
  * [18] A. D. Jones _et al._ , “IM3/HyperFACETS Thermodynamic Global Warming (TGW) Simulation Datasets (v1.0.0) [Data set],” _MSD-LIVE Data Repository_ , 2022, https://doi.org/10.57931/1885756.
  * [19] A. D. Jones _et al._ , “Continental United States climate projections based on thermodynamic modification of historical weather,” submitted.
  * [20] J. G. Powers _et al._ , “The weather research and forecasting model: Overview, system efforts, and future directions,” _Bulletin of the American Meteorological Society_ , vol. 98, no. 8, pp. 1717–1737, Aug 2017.
  * [21] W. C. Skamarock _et al._ , “A description of the advanced research WRF model version 4,” Tech. Rep., 2019.
  * [22] C. R. McGrath _et al._ , “tell: a python package to model future total electricity loads in the united states,” _Journal of Open Source Software_ , vol. 7, no. 79, p. 4472, 2022. [Online]. Available: https://doi.org/10.21105/joss.04472
  * [23] Electric Vehicle Infrastructure Projection Tool (EVI-Pro) Lite. Accessed: May 15, 2023. [Online]. Available: https://tinyurl.com/bdfjjfuu
  * [24] Fleet DNA: Commercial Fleet Vehicle Operating Data. Accessed: Jun 16, 2022. [Online]. Available: https://tinyurl.com/3u299eek
  * [25] Passenger boarding and all-cargo data for U.S. airports. Accessed: Feb 14, 2023\. [Online]. Available: https://tinyurl.com/yeykhjtw
  * [26] County transportation profiles. Accessed: Feb 14, 2023. [Online]. Available: https://www.bts.gov/ctp
  * [27] U.S. Single family homes - statistics & facts. Accessed: May 17, 2023. [Online]. Available: https://tinyurl.com/4ecrcepr
  * [28] D. Smith _et al._ , “Medium-and heavy-duty vehicle electrification: An assessment of technology and knowledge gaps,” Oak Ridge National Lab.(ORNL), Oak Ridge, TN (United States), Tech. Rep., 2020.


# Transfer of energy from flexural to torsional modes for the Fish-bone
Suspension bridge

Clelia Marchionna  Dipartimento di Matematica del Politecnico, Piazza Leonardo
da Vinci 32, 20133 Milano, ITALY  and  Stefano Panizzi  Dipartimento di
Scienze Matematiche, Fisiche e Informatiche, Parco Area delle Scienze 53/A,
43126 Parma, ITALY

###### Abstract.

We consider a conservative coupled oscillators systems which arises as a
simplified model of the interaction of flexural and torsional modes of
vibration along the deck of the so-called fish-bone [4] model of suspension
bridges. The elastic response of the cables is supposed to be asymptotically
linear under traction, and asymptotically constant when compressed (a
generalization of the slackening regime). We show that for vibrations of
sufficiently large amplitude, transfer of energy from flexural modes to
torsional modes may occur provided a certain condition on the parameters is
satisfied. The main result is a non-trivial extension of a theorem in [15] to
the case when the frequencies of the normal modes are no more supposed to be
the same. Several numerical computations of instability diagrams for various
slackening models respecting our assumptions are presented.

Keywords: Suspension bridges, torsional instability, Hill equation,
instability tongues, coupled oscillators

Mathematics Subject Classification: Primary: 74K30; Secondary: 34C15, 34C25,
34B30, 35G61, 35Q74

##  1\. Introduction

In this paper, which is a completion of a previous work [15], we face the
question of energy exchange between torsional and longitudinal modes in a
generalized energy conserving suspension bridge model, the so-called _fish-
bone_ model, proposed by K.S. MOORE [21], in which the elastic response of the
cables is asymptotically linear under traction, and asymptotically constant
when compressed (a generalization of the slackening regime).

We follow the line of research carried out by F. GAZZOLA and coworkers in a
series of papers [2, 4, 10] according to which internal nonlinear resonances
may occur even when the aeroelastic coupling is disregarded. To analyze the
onset of resonances, the PDEs system (infinite dimensional system) is reduced
to a coupled system of nonlinear ODEs, by projecting the infinite dimensional
phase space on a two dimensional subspace through an approximate Galerkin
technique. The ODEs system is then a simplified model of the interaction of
flexural and torsional modes of vibration along the deck of the bridge.

The same problem was addressed in a previous paper [15], in which only the
interaction between the modes corresponding to the fundamental flexural and
torsional frequencies was considered. In the present work, we extend the study
to the more general case in which modes of different frequencies can interact.

The main result is Theorem 1 in which we prove that for deflections of
sufficiently large amplitude, pure flexural and periodic vibrations are
unstable, provided a certain condition on the structural parameters of the
bridge is satisfied, so that transfer of energy from flexural modes to
torsional modes may occur even when no external forces are applied.

Our technique follows the lines of previous literature: we consider the
differential of the Poincaré map relative to the torsional component of the
system, around a periodic purely flexural solution. As is well-known this
leads to a Hill equation; then, as the energy of the system tends to infinity,
we compute the limit equation, which turns out to be simple enough to allow
the analysis of its stability via the classical Floquet theory.

If the transverse and torsional modes are labeled by $j$ and $k$ respectively,
an unexpected property holds when $j$ is even: the instability condition at
large deflections of Theorem 1 is never satisfied, as the limit system is an
uncoupled system of linear, constant coefficients oscillators. Consequently,
every solution of the limit system is periodic thus stable, so that purely
flexural solutions are asymptotically linearly stable.

The extension of the result in [15] is non-trivial because 2 main steps are
involved: the quite subtle derivation of the $C^{1}$-regularity of the
projected system under the relaxed regularity assumption ($\bf S_{0}$) on the
slackening function; the computation of the limit Hill equation at large
energies, which requires several technicalities.

In order to obtain a more general picture of the stability of the flexural
component, in Section 3 we look at the Hill equation relative to the torsional
component as dependent on two parameters. The first natural parameter is the
maximum elongation of the purely flexural solution, which in Theorem 1 tends
to $\infty$; we note that the dependence on such parameter is nonlinear. The
choice of the second parameter is motivated by mathematical reasons, as the
equation naturally presents a spectral parameter which, from the structural
point of view, essentially corresponds to the torsional behavior of the
bridge. In this way, we can compare the results related to our problem with
the instability diagrams in the literature for classical two-parameters Hill
equations. In particular, we consider a slackening model that has minimal
regularity with respect to our requests, using both academic and real
parameters corresponding to Tacoma narrow bridge; we numerically draw the
corresponding instability diagrams, pointing out which properties are
preserved and which diverge from the classical ones; some interesting
mathematical aspects are highlighted, such as the presence of resonance
pockets typical of some two parameters Hill equations, e.g. the multi-step
Meissner equation.

The remainder of this section is devoted to the presentation of the PDEs
model, and its reduction to the ODEs system. In Section 2 we provide the proof
of the main instability result. The work ends with three appendices: in
Appendix A we prove the regularity of the ODEs system; the computation of the
limit system at large energies is presented in Appendix B; finally, Appendix C
provides the explicit computation of the instability discriminant for a
piecewise linear system leading to a multi-step Meissner equation.

###  1.1. The bridge model

We briefly recall the PDEs system. The dynamics of the midline of the deck,
modeled as an Euler-Bernoulli beam of length $L$ and width $2l$, is coupled
with the elastic response of the suspension cables acting on the side ends of
the deck. The cross section of the deck is assumed to be a rigid rod with mass
density $\rho$, length $2l$ and thickness small with respect to $l$; $Y(x,t)$
is the vertical downward deflection of the midline of the deck with respect to
the unloaded state, $\Theta(x,t)$ is the angle of rotation of the deck with
respect to the horizontal position. The corresponding PDEs system is given by,

(1) $\begin{cases}Y_{tt}+\frac{EI}{\rho
S}\,Y_{xxxx}+f(Y+l\sin\Theta)+f(Y-l\sin\Theta)=0,\\\\[4.30554pt]
\Theta_{tt}-\frac{GK}{\rho
J}\Theta_{xx}+\,\frac{Sl}{J}\cos\Theta\left[f(Y+l\sin\Theta)-f(Y-l\sin\Theta)\right]=0,\end{cases}$

complemented with hinged boundary conditions:

(2) $Y(0,t)=Y(L,t)=Y_{xx}(0,t)=Y_{xx}(L,t)=0,\quad\Theta(0,t)=\Theta(L,t)=0.$

The other constant parameters are: $S$ the cross section area, $I$ the planar
second moment of area with respect to the plane $Y=0$, $J$ the polar second
moment of area with respect to the $x$-axis and $E$ and $G$ respectively the
Young modulus and the shear modulus, K the torsional constant.

The restoring force $f$ exerted by the hangers is applied to both extremities
of the deck whose displacements from the unloaded state are given by $Y\pm
l\sin\Theta$. No external forces, except gravity, are taken in account.

In the classical slackening regime, the hangers behave as linear springs of
elastic constant ${\rm k}>0$ if stretched and do not exert restoring force if
compressed. A first model in which the system (1) acts as a linear, non
coupled system for small displacements, was proposed by K.S MOORE and P.J.
McKENNA [18, 21] (hereforth MMK model) assuming for $f$ the following
expression ($g$ is the gravity),

(3) $f(r)=m\left[(r+r_{0})^{+}-r_{0}\right],\quad r_{0}=\rho Sg/2{\rm k},\quad
m=\frac{{\rm k}}{\rho S}.$

Subsequently many other forms for $f$ have been proposed in [4, 15, 16, 19],
some of these are nonlinear and smooth in a neighborhood of the origin, making
instability feasible even at low energies. Two significant examples are ($h$
is a positive constant):

(4) $f(r)=mr+\sqrt{(mr)^{2}+h^{2}}\,\,-h,$ (5) $f(r)=h(e^{mr/h}-1),$

Throughout the paper we assume that the function $f$ satisfies the following
mild regularity condition:

_Assumption_ ($\bf S_{0}$)

  * a)

$f$_is continuous, strictly increasing, and $f(0)=0$_;

  * b)

$f$ _is piecewise $C^{1}$_, that is _its derivative is continuous with the
exception of a finite (eventually empty) set of points $r_{1}<r_{2}<...<r_{n}$
not including zero in which there exist the finite limits: _

$\lim_{r\rightarrow r_{i}^{\pm}}f^{\prime}(r);$

  * c)

$m:=f^{\prime}(0)>0$.

In most cases the elastic response of the cables is supposed to be
asymptotically linear under traction, and asymptotically constant when
compressed, thus it is natural to assume at least one of the following
conditions:

  * ($\bf S_{1}$)

$\lim_{r\rightarrow-\infty}f^{\prime}(r)=0$

  * ($\bf S_{2}$)

$M:=\lim_{x\rightarrow+\infty}f^{\prime}(r)>0.$

We note that both (3) and (4) satisfy $(\bf S_{0})$ \- $(\bf S_{1})$ \- $(\bf
S_{2})$, unlike the function (5) which does not satisfies condition $(\bf
S_{2})$.

The problem (1)-(2) is well-posed in appropriate Sobolev spaces [4, 8], and
enjoys two properties mostly relevant for our purposes: the total energy is
conserved over time; it admits _pure flexural solutions_ , that is motions in
which the cross sections of the deck remain horizontal at all times so that no
torsional vibrations occur.

###  1.2. The ODEs system

For small torsional angles, it is very convenient to replace the system (1)
with a pre-linearized one, see [4, 15]. By the usual approximation:
$\sin\Theta\sim\Theta$, $\cos\Theta\sim 1$, and by setting $Z=l\Theta$, the
system (1) reduces to

(6) $\begin{cases}\,Y_{tt}+\frac{EI}{\rho
S}\,Y_{xxxx}+f(Y+Z)+f(Y-Z)=0\\\\[4.30554pt] \,Z_{tt}-\frac{GK}{\rho
J}\,Z_{xx}+\,\frac{l^{2}S}{J}\left[f(Y+Z)-f(Y-Z)\right]=0,\end{cases}$

and as far as the scope of this paper is concerned, nothing changes starting
from the system (1) or from (6).

Our _ansatz_ is that, after a suitable rescaling of the space variable, the
displacements can be reasonably well approximated by the $j-k$ mode of
vibration, that is,

$Y(x,t)\simeq y_{j}(t)\sin(jx),\quad Z(x,t)\simeq z_{k}(t)\sin(kx),\quad 0\leq
x\leq\pi.$

Then, through a Galerkin projection, the PDEs system reduces to a coupled
oscillators system which, dropping the indexes $j-k$, reads as follows:

(7) $\begin{cases}\ddot{y}+\alpha j^{4}y+\psi_{1}(y,z)=0\\\\[4.30554pt]
\ddot{z}+\beta k^{2}z+\gamma\psi_{2}(y,z)=0,\end{cases}$

with structural parameters, $\alpha=EI\pi^{4}/\rho SL^{4}$,
$\beta=GK\pi^{2}/\rho JL^{2}$, $\gamma=l^{2}S/J$, and nonlinear coupling
terms,

(8)
$\displaystyle\psi_{1}(y,z)=\frac{2}{\pi}\int_{0}^{\pi}\left[f(y\sin(jx)+z\sin(kx))+f(y\sin(jx)-z\sin(kx))\right]\sin(jx)\,dx,$
(9)
$\displaystyle\psi_{2}(y,z)=\frac{2}{\pi}\int_{0}^{\pi}\,\left[f(y\sin(jx)+z\sin(kx))-f(y\sin(jx)-z\sin(kx))\right]\sin(kx)\,dx.$

If we define

$\Psi(y,z)\,=\,\frac{2}{\pi}\int_{0}^{\pi}\left[\mathcal{F}(y\sin(jx)+z\sin(kx))+\mathcal{F}(y\sin(jx)-z\sin(kx))\right]\,dx,$

where $\mathcal{F}(r)=\int_{0}^{r}f(s)\,ds$, we have
$\psi_{1}\,=\,\partial\Psi/\partial y$, $\psi_{2}\,=\,\partial\Psi/\partial
z$, so that the system (7) admits a conserved energy,

(10)
$\mathcal{E}(y,\dot{y},z,\dot{z})\,=\,\frac{\dot{y}^{2}}{2}+\frac{\dot{z}^{2}}{2\gamma}+\frac{\alpha
j^{4}}{2}y^{2}+\frac{\beta k^{2}}{2\gamma}z^{2}+\Psi(y,z).$

Note that, under Assumption ($\bf S_{0}$), $\Psi(y,z)$ is nonnegative. As a
consequence all solutions of (7) are global and bounded.

Since $\psi_{2}(y,0)\equiv 0$, the system (7) admits periodic pure flexural
solutions, that is solutions of the form $y=u(t)$, $z\equiv 0$ with $u(t)$
periodic. We consider such solutions as parametrized by the initial
displacement, and we define $u=u(t;q)$ as the solution of the initial value
problem,

(11) $\ddot{u}+\alpha j^{4}u+2f_{j}(u)=0\qquad u(0)=q,\quad\dot{u}(0)=0.$

where

(12) $f_{j}(r):=\frac{1}{2}\,\psi_{1}(r,0)=\frac{2}{\pi}\int_{0}^{\pi}f(r\sin
jx)\sin jx\,dx.$

Assuming for the moment that $f\in C^{1}$, the linearization at a fixed energy
level (iso-energetic linearization) of the system around the periodic orbit
$(u(\cdot,q),0)$ yields, for the torsional component, the Hill equation (see
e.g. [6, 15] for details),

(13) $\ddot{v}+\left(\beta k^{2}+2\gamma g_{j,k}(u(t;q))\right)v=0,$

in which we have set,

(14) $g_{j,k}(r)=\frac{1}{2}\,\frac{\partial\psi_{2}}{\partial
z}(r,0)=\frac{2}{\pi}\int_{0}^{\pi}f^{\prime}(r\sin jx)\sin^{2}kx\ dx.$

The problem we want to address is the stability of solutions of the Hill
equation (13). It is worth noting that in the case $j=k$, we have
$g_{j,j}(r)=f^{\prime}_{j}(r)$, and the linearized system (11)-(13) is the
same as the one studied in [15]. This is no longer true if $j\neq k$. In [15],
under the assumptions $\bf(S_{0})$ \- $\bf(S_{1})$ \- $\bf(S_{2})$, we
established a condition depending on a set of 3 parameters under which the
flexural motions are unstable provided the energy parameter $q$ is
sufficiently large. The next section is devoted to the main result of the
present paper which is a non-trivial extension of the result in [15] to the
case $j\neq k$.

##  2\. Instability of pure flexural $j-k$ modes at high energies

The stability analysis of (13) is carried out by means of Floquet’s theorem,
see [7, 14]. We recall the definition of the stability discriminant
$\Delta=\Delta(q)$ of the Hill equation. Let $\tau(q)$ be the period of the
solution $u$ of the problem (11), and let $v_{0}(t)$, $v_{1}(t)$ be the
solutions of (13) corresponding to initial data

$v_{0}(0)=1,\,\dot{v}_{0}(0)=0,\quad v_{1}(0)=0,\,\dot{v}_{1}(0)=1;$

then

$\Delta(q)=v_{0}(\tau(q))+v^{\prime}_{1}(\tau(q)).$

If $|\Delta|>2$ the non-trivial solutions of the Hill equation are unbounded,
if $|\Delta|<2$ are all bounded. In the case when $\Delta=2$ there exists at
least a non trivial $\tau(q)$-periodic solution, when $\Delta=-2$ there exists
at least a non-trivial $2\tau(q)$-periodic solution.

The main result of this section consists in the the computation of

$\Delta_{\infty}:=\lim_{q\rightarrow\infty}\Delta(q),$

in the case when the elastic response of the cables is asymptotically linear,
i.e under assumptions $\bf(S_{1})$ \- $\bf(S_{2})$. Referring to the
characterization of instability $|\Delta_{\infty}|>2$, we can establish a
condition, depending on $j$, $k$, and the structural parameters, for which
there is instability for sufficiently high energies. The proof mimics that of
Theorem 1 in [15]. First we are able to compute the limit system of (11)-(13)
as $q$ goes to $+\infty$; it turns out that the Hill equation of the limit
system is a two-step potential (Meissner equation) that can be integrated
explicitly; the condition (15) expresses the instability condition for $q$
sufficiently large through the discriminant $\Delta_{\infty}$ of the limit
equation.

######  Theorem 1.

Assume that the function $f$ satisfies the conditions $\bf(S_{0})$ \-
$\bf(S_{1})$ \- $\bf(S_{2})$.

Assume that $\,j\,$ is odd, and let the constants $\omega_{\pm}$, $A_{\pm}$,
$\phi_{\pm}$, be defined as follows:

$\omega_{\pm}^{2}=\alpha j^{4}+M\,\left(1\pm\frac{1}{j}\right),$
$A_{\pm}^{2}=\beta k^{2}+\gamma
M(1\pm\epsilon_{j,k}),\quad\epsilon_{j,k}=\frac{1}{j}-\frac{\tan(\pi
k/j)}{k\pi},$ $\phi_{\pm}=\frac{A_{\pm}}{\omega_{\pm}}\,\pi,\qquad
a=\frac{A_{+}}{A_{-}}.$

Then, if the following condition holds true,

(15)
$\frac{|\Delta_{\infty}|}{2}:=\left|\cos\phi_{+}\,\cos\phi_{-}-\frac{a+a^{-1}}{2}\,\sin\phi_{+}\,\sin\phi_{-}\,\right|>1.$

there exists $\,q_{0}\,$ such that, if $\,q>q_{0}$, the pure flexural periodic
solution $(u(t;q),0)$ of the (non-linear) system (7) is unstable.

From the above expression of $A_{\pm}^{2}$ it is not obvious that they are
positive. Actually the two quantities $1\pm\epsilon_{j,k}$ result as the two
possible values ($+$ for $r>0$, and $-$ for $r<0$) of the positive integral
$\frac{4}{\pi}\int_{0}^{\pi}H(r\sin jx)\sin^{2}kx\ dx$, see (18), (20) below,
and Appendix B.

######  Remark 1.

In Theorem 1 we have considered only the case of odd $j$ as for even $j$ the
condition (15) is never satisfied. Indeed, thanks to (19) the limit system is
simply a decoupled linear system that we can write explicitly as follows,

$\begin{cases}\ddot{U}_{\infty}+(\alpha j^{4}+M)U_{\infty}=0\\\\[4.30554pt]
\ddot{v}_{\infty}+(\beta k^{2}+\gamma M)v_{\infty}=0.\end{cases}$

###  2.1.  Linearized $j-k$ mode: technical tools

First of all we need a regularity result which is crucial for the
linearization process of the system (7) around a pure flexural solution.

######  Proposition 2.1.

If the function $f$ satisfies assumption ($\bf{S_{0}}$), then the functions
$\psi_{1}$, $\psi_{2}$, as defined in (8)-(9), are of class
$C^{1}(\mathbb{R}^{2})$.

The proof of this Proposition is given in Appendix A.

Then we list some general facts about the functions $f_{j}(r)$, $g_{j,k}(r)$
defined in (12),(14). Given any function $\,f\,$, accordingly with [15], we
define the integral transform,

(16) $\tilde{f}(r):=\frac{2}{\pi}\int_{0}^{\pi}f(r\sin x)\sin x\,dx.$

We denote by $\,f_{e}(x)=\frac{1}{2}\left(f(x)+f(-x)\right)$,
$\,f_{o}(x)=\frac{1}{2}\left(f(x)-f(-x)\right)$ even and odd parts of a
function respectively. We have the following,

######  Lemma 1.

If the function $\,f\,$ satisfies the first two conditions of the assumption
($\bf S_{0}$), then $f_{j}\in C^{1}(\mathbb{R})$, $f_{j}(0)=0$ and
$f^{\prime}_{j}\geq 0$.

Moreover we have,

(17) $f_{j}=\tilde{f_{o}}\quad(\text{even }j),\qquad
f_{j}=\tilde{f_{o}}+\frac{1}{j}\,\tilde{f_{e}}\quad(\text{odd }j).$

###### Proof.

The regularity and the sign of $\,f^{\prime}_{j}\,$ follow from the analogous
properties of $\tilde{f}$ proved in [15], condition $\mathbf{(\tilde{H})}$. To
prove (17), we use the elementary identity $\,\sin(z+h\pi)=(-1)^{h}\sin z$, to
compute

$\displaystyle\int_{0}^{\pi}f(r\sin jx)\sin jx\,dx$ $\displaystyle=$
$\displaystyle\frac{1}{j}\int_{0}^{j\pi}f(r\sin z)\sin z\,dz$ $\displaystyle=$
$\displaystyle\frac{1}{j}\sum_{h=0}^{j-1}(-1)^{h}\int_{0}^{\pi}f((-1)^{h}r\sin
z)\sin z\,dz.$

By collecting the signs, and by distinguishing the two cases even or odd $j$,
we obtain (17). ∎

######  Lemma 2.

If $j$ is even, then $g_{j,k}(r)$ is an even function for every $k$. More
precisely, we have

$g_{j,k}(r)=\frac{2}{\pi}\int_{0}^{\pi}f_{o}^{\prime}(r\sin
jx)\sin^{2}kx\,dx.$

###### Proof.

If $j$ is even, $f^{\prime}(r\sin jx)(\sin kx)^{2}$ has period $\pi$, then

$g_{j,k}(r)=\frac{2}{\pi}\int_{0}^{\pi}f^{\prime}(r\sin
jx)\sin^{2}kx\,dx=\frac{2}{\pi}\int_{-\pi/2}^{\pi/2}f^{\prime}(r\sin
jx)\sin^{2}kx\,dx.$

It follows that

$g_{j,k}(r)=\frac{2}{\pi}\int_{-\pi/2}^{\pi/2}(f_{e}^{\prime}(r\sin
jx)\sin^{2}kx+f_{o}^{\prime}(r\sin jx)\sin^{2}kx\,)\,dx$

and the integral of the first term is null since $f_{e}^{\prime}$ is odd. ∎

We make an observation to highlight some differences that may exist between
even and odd harmonics of the flexural component. If the function $f$ is as in
example (4), its odd part is linear, precisely we have $f_{o}(r)=mr$. Then, if
$j$ is even, the equation (13) reduces to the linear oscillator,

$\ddot{v}+(\beta k^{2}+2\gamma m)v=0.$

Since $\beta k^{2}+2\gamma m>0$, the $j-k$ modes are always linearly stable.
This is no more true for odd $j$, see formula (17).

We end this section by introducing the functions which characterize the limit
system as $q\rightarrow\infty$ ($(\cdot)^{+}=$ positive part, $H(\cdot)$
Heaviside function). We set

(18) $h_{j}(r)=\frac{2M}{\pi}\int_{0}^{\pi}(r\sin jx)^{+}\sin jx\ dx,\qquad
s_{j,k}(r)=\frac{2M}{\pi}\int_{0}^{\pi}H(r\sin jx)\sin^{2}kx\ dx,$

then, by setting $r_{0}=0$ and $\theta=0$ in Lemma 5 in Appendix B, we obtain
their explicit expressions:

(19) $h_{j}(r)=\frac{M}{2}r,\quad s_{j,k}(r)=\frac{M}{2}\qquad(\text{even
}j);$ (20) $h_{j}(r)=\frac{M}{2}\left[1+\frac{{\rm sign}(r)}{j}\right]r,\quad
s_{j,k}(r)=\frac{M}{2}\left[1+{\rm
sign}(r)\left[\frac{1}{j}-\frac{\tan(k\pi/j)}{k\pi}\right]\right]\quad(\text{odd
}j).$

###  2.2. Proof of Theorem 1

###### Proof.

In [15] we proved the same facts for the simple mode with $j=k=1$, therefore
here we just outline the main points emphasizing a few differences.

Let $u$ be the solution of the problem (11), and $v$ be the solution equation
(13) with fixed initial data $v(0)=a,\,\dot{v}(0)=b$.

We rescale $u$ by setting $U_{q}(t)=u(t)/q$, then the system becomes

(21) $\begin{cases}\ddot{U}_{q}+\alpha
j^{4}U_{q}+2\displaystyle{\frac{f_{j}(qU_{q})}{q}}=0\\\\[4.30554pt]
\ddot{v}_{q}+(k^{2}\beta+2\gamma g_{j,k}(qU_{q}(t)))v_{q}=0\end{cases}$

with initial data $U_{q}(0)=1,\,\dot{U}_{q}(0)=0$,
$v_{q}(0)=a,\,\dot{v}_{q}(0)=b$.

Then we introduce the limit system of (21), as $q\rightarrow\infty$,

(22) $\begin{cases}\ddot{U}_{\infty}+\alpha
j^{4}U_{\infty}+2\displaystyle{h_{j}(U_{\infty})}=0\\\\[4.30554pt]
\ddot{v}_{\infty}+(\beta k^{2}+2\gamma
s_{j,k}(U_{\infty}(t))v_{\infty}=0\end{cases}$

with initial data $U_{\infty}(0)=1,\,\dot{U}_{\infty}(0)=0$,
$v_{\infty}(0)=a,\,\dot{v}_{\infty}(0)=b$.

We claim that, as $q\to+\infty$, the solutions of (21)-(22) satisfy the limit
relations:

  * i)

$U_{q}\to U_{\infty}\quad$ in $C^{0}(\mathbb{R})$;

  * ii)

$v_{q}\to v_{\infty}\quad$ in $C^{1}([0,T])$, for every $T>0$;

  * iii)

$\tau_{q}\to\tau_{\infty},\quad$ where $\tau_{q}$, $\tau_{\infty}$ are the
periods of $U_{q}$, $U_{\infty}$ respectively.

Proof of i): Let $\tilde{f}(r)$ be as in (16), then in [15] we proved that
$\lim_{q\to+\infty}\tilde{f}(qr)/q=\tilde{h}(r)$, uniformly on compact sets.
Thanks to (17), in the same way we obtain
$\lim_{q\to+\infty}f_{j}(qr)/q=h_{j}(r)$, uniformly on compact sets. Owing to
a classical continuous dependence theorem ([11], Th. 3, Ch.XV, p.297), we get
that $U_{q}$ converges uniformly on $\mathbb{R}$ to $U_{\infty}$.

Proof of ii): Let us consider the Hill equation in (22). By Lemma 5.1 in [15]
(see also [22],) we have to prove that $g_{j,k}(qU_{q}(t))\to
s_{j,k}(U_{\infty}(t))$ in $L^{1}[0,T]$.

In the case when $\sin jx\neq 0$, we have $U_{\infty}(t)\neq 0$, thus by the
uniform convergence of $U_{q}$, for $q$ sufficiently large, the sign of
$\,qU_{q}\sin jx\,$ is the same as that of $\,U_{\infty}\sin jx$. Thanks to
the assumptions (${\bf S_{1}}$)-(${\bf S_{2}}$), it follows that

$\lim_{q\rightarrow+\infty}f^{\prime}(qU_{q}(t)\sin jx)=MH(U_{\infty}(t)\sin
jx).$

Since $U_{\infty}(t)\neq 0$ almost everywhere in $[0,T]$, and $s_{j,k}$ is
bounded, we get the required convergence of $g_{j,k}(qU_{q})$ to
$s_{j,k}(U_{\infty})$.

The proof of iii) is the same as in Proposition 5.2 of [15].

Finally, we come to the computation of the limit period $\,\tau_{\infty}\,$
and discriminant $\,\Delta_{\infty}$. Recall that $j$ is odd, then the limit
system is a coupled (nonlinear) system with step coefficients $h_{j}(r)$,
$s_{j,k}(r)$ defined by formula (20). It turns out that the limit Hill
equation (22) is a two-step Meissner equation that can be integrated
explicitly. First of all we note that the function $U_{\infty}(t)$ is even, so
that we can fix our attention only on the half period. If we divide the
interval $[0,\tau/2]$ in two subintervals $I^{+}=[0,t_{0}]$, where
$U_{\infty}\geq 0$, and $I^{-}=[t_{0},t_{1}]$, where $U_{\infty}\leq 0$, we
easily obtain

$t_{0}=\frac{\pi}{2\omega_{+}},\qquad
t_{1}=\frac{\pi}{2\omega_{-}},\qquad\tau=2(t_{0}+t_{1}).$

It follows that the coefficient of the equation (22) is given by (here
$\mathbbm{1}_{I}(t)$ denotes the indicator function of $I$),

$\beta k^{2}+2\gamma
s_{j,k}(U_{\infty}(t))=A^{2}_{+}\mathbbm{1}_{I^{+}}(t)+A^{2}_{-}\mathbbm{1}_{I^{-}}(t)\qquad(0\leq
t\leq\tau/2).$

A straightforward computation (see e.g. [15] or Appendix C), shows that

$\frac{\Delta_{\infty}}{2}=v_{\infty}(\tau)=\cos\phi_{+}\,\cos\phi_{-}-\frac{a+a^{-1}}{2}\,\sin\phi_{+}\,\sin\phi_{-}.$

By the instability characterization of Hill equations, if
$|\Delta_{\infty}|>2$, then the pure flexural solutions $u(t;q)$ are linearly
unstable (thus unstable) for sufficiently large $q$. ∎

##  3\. The MMK model: instability tongues and some numerical results

A classical problem in the theory of Hill equations consists in studying the
parametric resonance of equations of the type

(23) $\ddot{v}+(\beta+qp(t)))v=0.$

Here $p(t)$ is a periodic function of period $\pi$ (just for fixing one), and
$\beta$, $q$ are real parameters. The general picture is well known and can be
briefly described as follows [14, ch. II, Th. 2.1] [7, ch. 2, Th. 2.3.1]. In
the $(q,\beta)$-plane a sequence of separated regions of instability
(instability tongues, Arnold’s tongues) emerge from the points $(0,n^{2})$,
$n=1,2,...$ on the $\beta-$axis which are bounded by two curves
$\beta=\beta_{n}^{\pm}(q)$. For values $(q,\beta)$ in the interior of the
tongues, i.e for $\beta^{-}_{n}(q)<\beta<\beta^{+}_{n}(q)$, all solutions of
(23) are unstable, while outside are stable. At the boundaries,
$\beta=\beta_{n}^{\pm}(q)$, we have at least a non-trivial $\pi$-periodic or
$\pi$-anti-periodic solution. For some $p(t)$, and for some exceptional values
of $q$, the curves $\beta^{\pm}(q)$ may intersect (co-existence of periodic
solutions) and form the so-called resonance pockets, as in the case of the
square wave case $p(t)=\text{sign }\cos(2t)$.

The classical, and most studied, example is the Mathieu equation,
$p(t)=\cos(2t)$, whose instability diagram was first drawn by B. VAN DER POL &
M. J. O. STRUTT in [25] and can be found on many books [20, 3]. In this case
we know the order of the instability tongues as $q\rightarrow 0$, the
asymptotic behavior of their width as $n\rightarrow\infty$, and many other
features, such as the fact that the width of the stability bands shrinks to
zero as $q\rightarrow\infty$, so that the areas of instability invade the
whole plane, see [13, 26].

Another case extensively studied is when $p(t)$ is a two-step function in the
interval $[0,\pi]$ (Meissner equation). The papers [12, 24] provide very
detailed results on the existence of resonance pockets and on the asymptotical
behavior of the stability boundaries.

In [17] we studied a generalization of the Mathieu equation,

(24) $\ddot{v}+(\beta+g(u(t;q)))v=0,$

in which the periodic coefficient $g(u(t;q))$ depends on the solution
$u=u(t;q)$ of an initial-value problem for a conservative second order
differential equation,

(25) $\ddot{u}+4u+f(u)=0,\qquad u(0)=q,\quad\dot{u}(0)=0,$

where $f(r)=O(r^{2})$, $r\rightarrow 0$. Note that the Mathieu equation
corresponds to the case $f=0$ in (25), and $g(u)=u$ in (24). The equation (24)
is again a Hill equation with two parameters that shares some aspects with the
Mathieu equation. Indeed the main result in [17] concerns the order of
tangency at $q=0$ of the instability tongues:

(26) $\beta_{n}^{+}(q)-\beta_{n}^{-}(q)=O(q^{n}),\quad q\rightarrow 0,$

which is the same as that of the Mathieu equation, at least if $f$ and $g$ are
real analytic functions near the origin. However, two features make its
analysis different (and more difficult): its period depends on the parameter
$q$, as it is half the period of $u(t;q)$ if $g$ is an even function, and is
the same otherwise; the dependence on $q$ in (24) is no more linear.

The system (11)-(13) has the same structure of (25)-(24), then for regular
$f$, $g$ the asymptotic formula (26) applies. Having at our disposal a result
for high energies (Theorem 1) and another one for low energies (formula (26)),
it would be interesting to have a general picture at least in the significant
case of the MMK restoring force (3), which is a simple non-trivial model
satisfying the Assumptions ($\bf S_{0}$)-($\bf S_{1}$)-($\bf S_{2}$), being
not smooth in a single point $r_{0}$.

We now come to the presentation of the models adopted in our numerical
simulations. We focus only on the system (11)-(13) in the case when function
$f$ is the MMK function (3), and on the system (29)-(30) below, which is a non
smooth variation of it in the case $j=k$. We performed our simulations both
with a simple set of parameters with the same order of magnitude (Figure 1),
and with the set as in [8, 9] where the mechanical parameters of the Tacoma
narrow bridge (TNB) and of other significant bridge models are listed,
including the MMK model, and several numerical experiments are performed.

The analytic expression of the functions $f_{j}$, owing to (17) in Lemma 1, is
known once we provide the projected form $\tilde{f}$ of $f$ in (3). This was
already done in our previous paper [15], and its closed form is given by,

(27) $\tilde{f}(r)=mr,\quad r\geq-
r_{0},\quad\tilde{f}(r)=-\frac{2}{\pi}mr_{0}\left[\frac{r}{r_{0}}{\rm
asin}\left(\frac{r_{0}}{r}\right)+\sqrt{1-\left(\frac{r_{0}}{r}\right)^{2}}\right],\quad
r\leq-r_{0}.$

As a consequence, for $j=k$, also the the functions $g_{j,j}=f^{\prime}_{j}$
(see Lemma 1) have a simple closed form. Much more complicated is the
derivation of the functions $g_{j,k}$ in (14) when $j\neq k$. Their complete
computation is provided in Appendix B.

In order to compare different functions, modeling the same slackening regime,
we introduce only for the case $j=k$ a non smooth variation of $f_{j}$, and
$f^{\prime}_{j}$ which maintains the same shape as the original MMK model but
which is adjusted in a way to have the same asymptotic behavior of $f_{j}$,
that is $\lim_{r\rightarrow-\infty}\tilde{f}(r)=-\frac{4}{\pi}mr_{0}$. More
precisely, we replace the function $\tilde{f}$ in (27) with the function

(28) $\bar{f}(r)=m\left[(r+4r_{0}/\pi)^{+}-4r_{0}/\pi\right]$

Then $f_{j}$ can be approximated, according to Lemma 1, with

$\bar{f}_{j}(r)=\frac{1}{2}\,\bar{f}_{o}(r),\qquad(\text{even
}j),\qquad\bar{f}_{j}(r)=\frac{1}{2}(\bar{f}_{o}(r)+\frac{1}{j}\bar{f}_{e}(r)),\qquad(\text{odd
}j),$

and we introduce the system,

(29) $\ddot{u}+\alpha j^{4}u+2\bar{f}_{j}(u)=0,\,\,u(0)=q,\,\,\dot{u}(0)=0$
(30) $\ddot{v}+(\beta j^{2}+2\gamma\bar{f}^{\prime}_{j}(u))v=0,$

Both $u$ and $v$ can be calculated explicitly and in Appendix C we provide the
formula for the instability discriminant $\Delta$ for every $(q,\beta)$. The
expression of $\Delta(q,\beta)$ is very complicated and is actually of little
help from the analytical point of view, but it can be used to represent very
quickly, using Matlab, the instability diagram in order to make a comparison
with the numerical results regarding the system (11)-(13).

Before proceeding with the discussion of the numerical results, let us fix the
starting points $\,\beta_{N}^{+}(0)=\beta_{N}^{-}(0)$ of the instability
tongues for our equations. Recall that we use $\beta$ as spectral parameter.

######  Proposition 3.1.

Let us suppose that the function $f$ satisfies the assumption ($\bf S_{0}$).
Then instability tongues for the Hill equation (13) stem from the values of
$\beta_{N}(0)$, $N=1,2...$, in the $\beta$-axis, given by

$\beta_{N}(0)=\frac{\alpha j^{4}+2m}{4k^{2}}\,N^{2}-\frac{2\gamma m}{k^{2}}.$

If $j$ is even the tongues corresponding to an odd index $N$ disappear, since
the actual period is half the period of $u(t;q)$.

We skip the proof which, after a translation of the spectral parameter
$\beta$, follows from the fact that, if $\tau(q)$ is the period of $u(t;q)$,
then its limit as $q\rightarrow 0$, is given by $\tau_{0}=2\pi/\sqrt{\alpha
j^{4}+2m}$. The different behavior for even or odd $j$ arises from the fact
that for even $j$, $g_{j,k}(u(t))$ is an even function and the minimal period
of the equation (13) is half the period $\tau(q)$ of $u(t;q)$.

Figure 1. Instability diagrams of systems (29)-(30), first line, and (11, 13),
second line. The fixed parameters are: $\alpha=1$, $\gamma=3$, $m=3$,
$r_{0}=1/3$. $\Delta(q,\beta)>2$ in t he light grey zones, and
$\Delta(q,\beta)<-2$ in the dark grey zones. Note that on the right column the
tongues corresponding to odd index vanish.

Figure 2. Instability diagrams of system (11)-(13) for the MMK function, and
$j\neq k$. The fixed parameters (TNB) are: $\alpha=8.0353\cdot 10^{-4}$,
$m=185.1$, $r_{0}=0.0265$. The actual value $\beta=8.1833\cdot 10^{-5}$ of the
TNB is highlighted.

In the first line of Figure 1 we show the diagram of the instability tongues
for the ”irregular approximation” (29)-(30) when $j=1$ (left), $j=2$ (right).
In the second line of Figure 1 the diagrams of the corresponding instability
tongues for the system (11)-(13), ($\tilde{f}$ as in (27)), with $j=k$ and
with the same parameters, in order to compare the diagrams with the first
line. It is worth noting that both systems are linear decoupled systems for
low values of $q$, then we have a trivial case of coexistence of periodic
orbits, since all solutions are periodic.

Coming to the specific features of the system (29)-(30), the true instability
tongues split from the points $(\frac{4}{\pi}r_{0},\beta_{N}(0))$, where
$\beta_{N}(0)$ is as in Proposition (3.1). We note the persistence of areas of
stability even for large values of $q$, so that their widths do not shrink to
$0$, as $q\rightarrow\infty$. This is more evident for greater values of $j$
and is in accordance with Theorem 1. In this regard, it is worth comparing the
instability diagram in Figure 1, and the one when $p(t)$ is a step function in
problem (23), see [24]. In both cases $\beta$ has the same meaning, but in the
first case the parameter $q$ acts both on the period and the width of the
steps of the Meissner equation, while in the second case, the parameters $q$
in (23) acts only on the height of the steps, as the period and width of the
steps are fixed. Then we note that the width at the splitting points
$(\frac{4}{\pi}r_{0},\beta_{N}(0))$, unlike (26), is $O(q)$, $q\rightarrow 0$.
This is due to the singularity of $\bar{f}_{j}$ in $r=4\pi/r_{0}$. Finally we
observe the onset of the so called resonance pockets, that are typical of some
Hill equations of Meissner type, see [12, 24].

The diagrams in the second line of Figure 1 are obtained by the numerical
solutions of the system (11)-(13). We used the Matlab solver ode23t with a
quite high reltol and absTol setting. Even though we used the closed formulas
for $f_{j}$, $g_{j,k}$ provided in Appendix B, a very long computation time
was needed to obtain an accurate value of $\Delta(q,\beta)$, and a reasonably
good quality diagram. Indeed, at every step in the grid value of $q$, a
precise evaluation of the period of $u(t;q)$ was necessary.

Figure 2 shows two instability diagrams of the system (11)-(13) respectively
for $j=3$, $k=2$ on the left, and $j=10$, $k=2$ on the right. In this case
$j\neq k$, so that the approximated system (29)-(30) is no more available.

More precisely, we decided to use the set of structural parameters of the
Tacoma Narrow bridge (TNB) reported in [8] that are derived from the Technical
Report (1941) of AMMANN et al. [1]. The mechanical parameters of the TNB are
characterized by very different orders of magnitude. Nonetheless the general
picture of the instability tongues (behavior at the splitting points and
resonance pockets) remains the same and seems to be a constant feature of the
MMK model.

The reports about the failure of the TNB show that the transfer from flexural
to torsional energy regarded most probably the $9-2$ or the $10-2$ modes (see
for example [10] for an extensive owerwiew). We checked all the modes $j$-2,
for $2\leq j\leq 10$, and $0<q\leq 1.6$ m (flexural oscillations of amplitude
of about one and a half meters were reported by witness, see for example again
[8], Appendix A). For our simplified model the tongues corresponding to an
even $N$ are relatively thick in the given interval for $q$, but stay far away
from the significant value of $\beta$. Whereas the tongues related to an odd
$N$ get very thin as $j$ grows, and also such modes do not present any
instability near the significant value of $\beta$. We point out that the fish-
bone model does not take into account the dynamics of the suspension cables
and their mechanical parameters, then our simulations are interesting only
from a mathematical point of view.

##  Appendix A Regularity of $\psi_{1}$ and $\psi_{2}$

######  Lemma 3.

Let $f$ be a real valued function such that $\,f\in C^{0}(\mathbb{R})$, $f\in
C^{1}(\mathbb{R}\setminus\\{r_{0}\\})$, and
$\quad\sup_{\mathbb{R}\setminus\\{r_{0}\\}}|f^{\prime}(x)|<\infty$,
$r_{0}\in\mathbb{R}$. Let $I$ be the closed interval $[a,b]$, $\,u\in
C^{1}(\mathbb{R}^{n+1})$, $w\in L^{\infty}(I)$, and define

$G(p)=\int_{I}f(u(x,p))w(x)\,dx,\qquad(p\in\mathbb{R}^{n}).$

If, for each $\,p\in\mathbb{R}^{n}$,

(31) $u(x,p)=r_{0},\text{ has (at most) a finite number of solutions}\quad
x\in I,$

then $\,G\in C^{1}(\mathbb{R}^{n})$, and the following differentiation formula
holds for every $p\in\mathbb{R}^{n}$,

$\bigtriangledown
G(p)=\int_{I}f^{\prime}(u(x,p))\,\bigtriangledown_{p}u(x,p)w(x)\,dx.$

###### Proof.

For simplicity we assume $n=1$. The general case $n\geq 1$ follows exactly in
the same way.

We set

$\Gamma=\\{(x,p)\in\mathbb{R}^{2}:\,u(x,p)=r_{0}\\},\qquad\Omega=\mathbb{R}^{2}\setminus\Gamma,$

and let $p_{0}\in\mathbb{R}$ be a number such that there exists a finite
number of points $(x_{j},p_{0})\in\Gamma$, $j=0,1,...,m$. Possibly by
decomposing the interval $I$ into a finite union of intervals $I_{j}$ in which
$u(x,p)=r_{0}$ has a single solution, we may assume that there exists a unique
point $(x_{0},p_{0})\in\Gamma$.

Let us fix $\varepsilon>0$ and define

$K_{\varepsilon}^{0}=\\{(x,p_{0}):\,x\in I,|x-x_{0}|\geq\varepsilon\\}.$

Since $K_{\varepsilon}^{0}$ is compact, its distance from the closed set
$\Gamma$ is strictly positive. Then there exists $\delta=\delta(\varepsilon)$,
such that

$K_{\varepsilon}^{\delta}=\\{(x,p):\,x\in
I,|x-x_{0}|\geq\varepsilon,|p-p_{0}|\leq\delta\\}\subseteq\Omega.$

By the continuity of $u$, we have either $u<r_{0}$ or $u>r_{0}$ on each
connected component of $K_{\varepsilon}^{\delta}$, thus $f\circ u\in
C^{1}(K_{\varepsilon}^{\delta})$.

Let us set

$\,g(p)=\int_{I}f^{\prime}(u(x,p))\,\frac{\partial u}{\partial
p}(x,p)w(x)\,dx,$

and split the function $G$ in the following way:

$G(p)=G_{1,\varepsilon}(p)+G_{2,\varepsilon}(p)\qquad(|p-p_{0}|\leq\delta),$

where

$G_{1,\varepsilon}(p)=\int_{x\in
I,|x-x_{0}|\geq\varepsilon}f(u(x,p))w(x)\,dx,\qquad
G_{2,\varepsilon}(p)=\int_{x\in I,|x-x_{0}|<\varepsilon}f(u(x,p))w(x)\,dx,$

so that $G_{1,\varepsilon}\in C^{1}([p_{0}-\delta,p_{0}+\delta])$ by the
standard rule for differentiation under the integral sign. As for the second
term we have,

$\displaystyle\left|\frac{G_{2,\varepsilon}(p_{0}+h)-G_{2,\varepsilon}(p_{0})}{h}-\int_{x\in
I,|x-x_{0}|<\varepsilon}f^{\prime}(u(x,p_{0}))\frac{\partial u}{\partial
p}(x,p_{0})w(x)dx\right|$ $\displaystyle\leq$ $\displaystyle\int_{x\in
I,|x-x_{0}|<\varepsilon}\left|\frac{f(u(x,p_{0}+h))-f(u(x,p_{0}))}{h}w(x)\right|+\left|f^{\prime}(u(x,p_{0}))\frac{\partial
u}{\partial p}(x,p_{0})w(x)\right|\,dx$

Since $f$ has bounded derivative (thus Lipschitz continuous), $w\in
L^{\infty}$, and $u\in C^{1}$, the sum of both terms in the last integral is
bounded above by some positive constant $M$. In conclusion, by considering the
incremental ratio of both $G_{1,\varepsilon}$, and $G_{2,\varepsilon}$ we get

$\limsup_{h\rightarrow
0}\left|\frac{G(p_{0}+h)-G(p_{0})}{h}-g(p_{0})\right|\leq 2M\varepsilon,$

for every $\varepsilon>0$, which proves that $G^{\prime}(p_{0})=g(p_{0})$.

Now we prove the continuity of $g(p)$. With the same notations as before, for
any fixed $\varepsilon>0$, we have

$g(p)=G^{\prime}_{1,\varepsilon}(p)+\int_{x\in
I,|x-x_{0}|<\varepsilon}f^{\prime}(u(x,p))\frac{\partial u}{\partial
p}(x,p)w(x)\,dx\qquad(|p-p_{0}|\leq\delta),$

where $G^{\prime}_{1,\varepsilon}\in C^{0}([p_{0}-\delta,p_{0}+\delta])$. As
for the other term, as before the boundedness of $f^{\prime}$, and the
regularity of $u$, yield

$\int_{x\in I,|x-x_{0}|<\varepsilon}\left|f^{\prime}(u(x,p))\frac{\partial
u}{\partial p}(x,p)w(x)\right|dx\leq
2M\varepsilon\qquad(|p-p_{0}|\leq\delta).$

Therefore we infer that

$|g(p)-g(p_{0})|\leq|G^{\prime}_{1,\varepsilon}(p)-G^{\prime}_{1,\varepsilon}(p_{0})|+4M\varepsilon,$

so that

$\limsup_{p\rightarrow p_{0}}|g(p)-g(p_{0})|\leq 4M\varepsilon,$

which concludes the proof of the Lemma.

∎

Proof of Proposition 2.1. We apply Lemma 3 to $f$ satisfying ($\bf{S_{0}}$),
$I=[0,\pi]$, $p=(y,z)\in\mathbb{R}^{2}$, $u(x,p)=y\sin(jx)\pm z\sin(kx)$,
$w(x)=\sin(jx)$ or $w(x)=\sin(kx)$. The lemma has a local character and it is
enough to verify the hypotheses in a neighborhood of each point
$p_{0}=(y_{0},z_{0})\in\mathbb{R}^{2}$, therefore we may assume that the
function $f$ is not differentiable at a single point $r_{0}$, and that has
bounded derivative elsewhere. The regularity of $u$, and $w$ being obvious, we
must verify the assumption (31) for $r_{0}\neq 0$ (recall that this was
required in ($\bf{S_{0}}$)). By contradiction, if for a fixed
$(y,z)\in\mathbb{R}^{2}$, the equation $u(x,y,z)=r_{0}$ were satisfied for an
infinite set of $x\in I$, since $x\mapsto u(x,y,z)$ then $u(x,y,z)=r_{0}$ for
every $x\in\mathbb{R}$. This follows by the unique continuation principle
applied to the analytic function $x\mapsto u(x,y,z)$. In particular we would
have $u(0,y,z)=0=r_{0}$, which is a contradiction.

We conclude this Appendix with a few remarks on the differentiation Lemma 3.

The lemma, in spite of its simplicity, cannot be derived from the traditional
Lebesgue theorem of differentiation under integral, as one can easily verify.
A simple example is provided by the function $\,G(p)=\int_{0}^{1}|x-p|\,dx$,
which is differentiable everywhere, and of class $C^{1}$. As a matter of fact,
Lemma 3 is more a regularity result than a sufficient condition to
differentiate under integral sign.

The condition (31) serves to our purposes but can be easily relaxed, for
example by assuming that the set $\\{x\in I:\,u(x,p)=r_{0}\\}$ has (at most) a
finite number of accumulation points for any fixed $p$. A trivial example in
which this last condition is violated, for $p=0$, is given by
$G(p)=\int_{0}^{1}|p|\,dx$, which of course is not differentiable at $p=0$. In
this regard, we point out the relevance in the application of the Lemma to the
functions $\psi_{1}$, and $\psi_{2}$, in which $u(x,p)=y\sin(jx)\pm
z\sin(kx)$, to the assumption $r_{0}\neq 0$. On the other hand, any relaxed
version of condition (31) is far from being necessary, as the simple example
$G(p)=\int_{0}^{1}|p^{3}|\,dx$ shows. In literature there are other more or
less classical differentiation results, see _e.g._ [23], but we have found
that for our purposes the verification of the hypotheses would have been more
difficult than the direct proof of the lemma.

##  Appendix B Explicit formulas for the MMK function

The non-linear terms of the system (7) and the periodic coefficient of the
Hill equation (13) are both defined as an integral. This fact becomes very
time consuming when it comes to drawing the instability tongues with Matlab,
as the solution of some ten thousands of systems is needed. So, in the case of
the MMK slackening model (3), we decided to provide the explicit computation
of $g_{j,k}(r)$ in (13), for every choice of $j$ and $k$. The result in Lemma
5 for $r_{0}=0$ is necessary for the computation of the limit system in
Section 2.

We need a simple computational lemma.

######  Lemma 4.

Let $j$, $k$ be two positive integers, then

(32)
$q(j,k):=\sum_{n=1}^{j}{\cos(\frac{2k}{j}n\pi)}\,=\,\begin{cases}j&\text{if
}\quad k/j\in\mathbb{N},\\\ 0&\text{otherwise}.\end{cases}$ (33)
$p(j,k):=\sum_{n=1}^{j}{(-1)^{n}\sin(\frac{2k}{j}n\pi)}=\begin{cases}0&\text{
if }j\text{ even },\\\ -\,\tan(\frac{k}{j}\pi)&\text{ otherwise}.\end{cases}.$

###### Proof.

Let us prove (32). In the case when $k/j\in\mathbb{N}$, we obviously have
$\,q(j,k)=j$; otherwise we get ($i=\sqrt{-1}$, $Re=$ real part),

$q(j,k)=\sum_{n=0}^{j-1}{\cos(\frac{2k}{j}n\pi)}=Re\,\sum_{n=0}^{j-1}e^{\frac{2k\pi
i}{j}n}=Re\,\frac{1-e^{2k\pi i}}{1-e^{2k\pi i/j}}=0.$

In order to prove (32), we observe that, if $\frac{2k}{j}\in\mathbb{N}$, then
$p(j,k)=0$. Otherwise, we have ($Im$ = imaginary part),

$p(j,k)=\sum_{n=0}^{j-1}{(-1)^{n}\sin(\frac{k}{j}2n\pi)}=Im\,\sum_{n=0}^{j-1}(-e^{\frac{k}{j}2\pi
i})^{n}=Im\,\frac{1-(-1)^{j}}{1+e^{\frac{k}{j}2\pi i}},$

thus $\,p(j,k)=0\,$ for $\,j$ even. When $j$ is odd, we get

$p(j,k)=-\,\frac{\sin(\frac{k}{j}2\pi)}{1+\cos(\frac{k}{j}2\pi)}=-\,\tan(\frac{k}{j}\pi).$

∎

The closed form of the function $\,g_{j,k}\,$, in the case of MMK function
(3), is given by $\,2m/\pi$ times the function $\,H_{j,k}(r)\,$ defined here
below.

######  Lemma 5.

Let $j$, $k$ be positive integers, $r_{0}\geq 0$, and let $H_{j,k}(r)$ be the
function defined as ($H[\cdot]=$ Heaviside step function),

$H_{j,k}(r):=\int_{0}^{\pi}H(r\sin jx+r_{0})\,\sin^{2}kx\
dx\qquad(r\in\mathbb{R}).$

Let us set $\;\theta(r)={\rm asin}(r_{0}/|r|)$, for $|r|\geq r_{0}$.

Then, for every $j$, $k$, we have

$\,H_{j,k}(r)=\frac{\pi}{2}\qquad\text{ for }\;|r|\leq r_{0},$

while for $|r|>r_{0}$ (note that for even $j$ the second row vanishes),

(34) $\displaystyle H_{j,k}(r)$ $\displaystyle=$
$\displaystyle\frac{\pi}{4}+\frac{\theta(r)}{2}-\frac{q(j,k)}{4k}\,\sin\left(\frac{2k\theta(r)}{j}\right)+$
$\displaystyle\,+\,\frac{1-(-1)^{j}}{2}\,\frac{r}{|r|}\left[\frac{\pi}{4j}-\frac{\theta(r)}{2j}+\frac{p(j,k)}{4k}\cos\left(\frac{2k\theta(r)}{j}\right)+\frac{1}{4k}\sin\left(\frac{2k\theta(r)}{j}\right)\right].$

###### Proof.

If $|r|\leq r_{0}$, we have $\,r\sin jx+r_{0}\geq 0$, thus
$H_{j,k}(r)=\int_{0}^{\pi}\sin^{2}kx\ dx=\pi/2$.

The other cases are much more involved. First of all, by changing the
integration variable, we write

$H_{j,k}(r)=\frac{1}{j}\int_{0}^{j\pi}H(r\sin z+r_{0})\sin^{2}(kz/j)\,dz;$

and observe that, for $\,z\in[0,2\pi]$, we have

$\,r\sin z+r_{0}\geq 0\;\text{ iff }\,z\in
B^{+}_{0}(r):=[0,\pi+\theta(r)]\cup[2\pi-\theta(r),2\pi]\qquad(r>r_{0}),$
$\,r\sin z+r_{0}\geq 0\;\text{ iff }\,z\in
B^{-}_{0}(r):=[0,\theta(r)]\cup[\pi-\theta(r),2\pi]\qquad(r<-r_{0});$

then we set $\,B^{\pm}_{n}(r)=B^{\pm}_{0}(r)+2n\pi$ (the translated sets).

We begin with $r>r_{0}$, and $j$ even. We get

$H_{j,k}(r)=\frac{1}{j}\sum_{n=0}^{j/2-1}\int_{B^{+}_{n}(r)}\sin^{2}(kz/j)\,dz;$

by direct computation of the integrals, we obtain

$\displaystyle\frac{1}{j}\sum_{n=0}^{j/2-1}\int_{B^{+}_{n}(r)}\sin^{2}(kz/j)\,dz$
$\displaystyle=$
$\displaystyle\frac{\pi}{4}+\frac{\theta}{2}+\frac{1}{4k}\sum_{n=0}^{j/2-1}\left(\sin(4kn\pi/j)-\sin(4k(n+1)\pi/j)\right)$
$\displaystyle+\frac{1}{4k}\sum_{n=0}^{j/2-1}(\sin(\frac{2k}{j}((2n+2)\pi-\theta))-\sin(\frac{2k}{j}((2n+1)\pi+\theta)).$

The first summation on the right hand side is telescopic and cancels out. The
last summation, by using the trigonometric addition formula, may be written as

$\frac{1}{4k}\left(\sum_{n=1}^{j}{(-1)^{n}\sin(\frac{2k}{j}n\pi)}\right)\cos(\frac{2k}{j}\theta)-\frac{1}{4k}\left(\sum_{n=1}^{j}{\cos(\frac{2k}{j}n\pi)}\right)\sin(\frac{2k}{j}\theta).$

Owing to (32), and (33), and putting together the various contributions, this
concludes the proof in the case $j$ even and $r>r_{0}$.

If $j$ is even and $r<-r_{0}$, we replace $B^{+}_{n}(r)$ with $B^{-}_{n}(r)$
and follow the same procedure, to obtain

$H_{j,k}(r)=\displaystyle{\frac{\pi}{4}+\frac{\theta}{2}-\frac{p(j,k)}{4k}\cos(\frac{2k}{j}\theta)-\frac{q(j,k)}{4k}\sin(\frac{2k}{j}\theta)}=\displaystyle{\frac{\pi}{4}+\frac{\theta}{2}-\frac{q(j,k)}{4k}\sin(\frac{2k}{j}\theta)}.$

Let us now consider the case when $j$ is odd. In the simple case $j=1$, the
direct computation of $H_{j,k}(r)$ is trivial, and yields

$H_{1,k}(r)=\frac{\pi}{2},\quad\text{ if }\,r>-r_{0};\qquad
H_{1,k}(r)=\theta(r)-\frac{\sin(2k\theta(r))}{2k}\quad\text{ if }\,r<-r_{0},$

which coincides with the formula (34) (although its verification is somehow
hidden).

For $j$ odd, and $j\geq 3$, and $r>r_{0}$, we get, again by following the same
procedure as above,

$\displaystyle H_{j,k}(r)$ $\displaystyle=$
$\displaystyle\frac{1}{j}\sum_{n=0}^{(j-3)/2}\int_{B^{+}_{n}(r)}\sin^{2}(\frac{k}{j}z)\,dz+\frac{1}{j}\int_{(j-1)\pi}^{j\pi}\sin^{2}(\frac{k}{j}z)\,dz$
$\displaystyle=\left(1+\frac{1}{j}\right)\frac{\pi}{4}+\left(1-\frac{1}{j}\right)\frac{\theta}{2}+\frac{p(j,k)}{4k}\cos(\frac{2k}{j}\theta)-\frac{1}{4k}(q(j,k)-1)\sin(\frac{2k}{j}\theta);$

while, for $r<-r_{0}$,

$\displaystyle H_{j,k}(r)$ $\displaystyle=$
$\displaystyle\frac{1}{j}\sum_{n=0}^{(j-1)/2}\int_{B^{-}_{n}(r)}\sin^{2}(\frac{k}{j}z)\,dz-\frac{1}{j}\int_{j\pi}^{(j+1)\pi}\sin^{2}(\frac{k}{j}z)\,dz$
$\displaystyle=\left(1-\frac{1}{j}\right)\frac{\pi}{4}+\left(1+\frac{1}{j}\right)\frac{\theta}{2}-\frac{p(j,k)}{4k}\cos(\frac{2k}{j}\theta)-\frac{1}{4k}(q(j,k)+1)\sin(\frac{2k}{j}\theta).$

Again owing to (32), and (33), this concludes the proof of the lemma. ∎

If $f(r)$ is the MMK function in (3), where $r_{0}>0$, then
$g_{j,k}(r)=\frac{2m}{\pi}H_{j,k}(r)$. We can verify by direct inspection in
(34) that this is a continuous function, as expected, because Proposition 2.1
says that $\frac{\partial\psi_{2}}{\partial z}(y,0)=2g_{j,k}(y)$ must be
continuous. If we set instead $r_{0}=0$, (34) gives us the explicit formula
for $s_{j,k}(r)=\frac{2M}{\pi}H_{j,k}(r)$, that is needed for defining the
limit system in Theorem 1. This function is not continuous in $r=0$ if $j$ is
odd.

This makes clear again that, if we weaken the condition b) in the assumption
($\bf S_{0}$), substituting the MMK function with the similar one
$f(r)=Mr^{+}$, the respective functions $\psi_{1}$, $\psi_{2}$ in (7) are no
more smooth.

##  Appendix C The discriminant for the approximation (28)

In this section we provide the computation of the instability discriminant
$\Delta$ for each fixed value of $\beta$ and $q$ for the Hill equation (30).
We find the explicit solution $u(t;q)$ of the equation (29), and its period
$\tau=\tau(q)$. The periodic coefficient in (30) is affected only by the slope
of the piecewise linear function $\bar{f}_{j}(r)$ whose values change at some
transition points of $u(t;q)$. It turns out that, at any fixed value of $q$,
(30) is a multi-step Hill equation or Meissner equation [5, 14].

We recall how to compute the discriminant of a Hill equation with a positive,
multi-step potential. Let the interval $[0,\tau]=\cup_{i=0}^{n}I_{i}$ be the
union of disjoint intervals, each having length $\,\Delta t_{i}$. Let us
consider a potential $Q(t)$ which is $\tau$-periodic positive and constant on
each subinterval $I_{i}$, that is,

$Q(t)=\sum_{i=0}^{n}A_{i}^{2}\,\mathbbm{1}_{I_{i}}(t)\qquad 0\leq t\leq\tau.$

The monodromy matrix, and the discriminant of the Hill equation
$\,\ddot{v}(t)+Q(t)v(t)=0$, are computed as follows, see [5] p. 12. If the
$L_{i}$’s are the transition matrices

$L_{i}=\begin{bmatrix}\cos(A_{i}\Delta t_{i})&\frac{1}{A_{i}}\sin(A_{i}\Delta
t_{i})\\\ -A_{i}\sin(A_{i}\Delta t_{i})&\cos(A_{i}\Delta
t_{i})\end{bmatrix},\quad i=0,1,...n,$

then we get

$M=L_{n}\,L_{n-1}...L_{1}\,L_{0},\qquad\Delta={\rm tr}M.$

To simplify some calculations, it is important to note that the discriminant
is invariant with respect to any cyclic permutation of the transition
matrices, so for instance $M=L_{2}\,L_{1}\,L_{0}$, and
$M^{\prime}=L_{1}\,L_{0}\,L_{2}$ are in general different matrices having the
same trace. In fact, they correspond to different translations in time of
$Q(t)$, which of course leave the discriminant unchanged.

In applying the previous formulas to the equation (30), the $A^{2}_{i}$’s
coefficients are known, being determined by the constant slopes of the
function $\bar{f}_{j}(r)$. It remains to compute the length of the intervals
$I_{i}$, and their ordering, modulo cyclic permutations.

First of all, we observe that, in both cases $j$ even or odd, we have

$\,\bar{f}_{j}(r)=mr,\quad\text{ if }\,|r|\leq\bar{r}:=\frac{4}{\pi}r_{0},$

therefore when the initial value $\,u(0)=q$ is less or equal to $\bar{r}$, the
equation (29) is linear, with solution $\,u(t;q)=q\cos(\omega t)$,
$\omega=\sqrt{\alpha j^{4}+2m}$. It follows that $\bar{f}_{j}^{\prime}$ equals
to $m$, and the Hill equation (30) reduces to a linear oscillator with
constant angular frequency $A=\sqrt{\beta j^{2}+2\gamma m}$. Thus it is stable
and $\Delta$ is simply the trace of the matrix $L$ with $\Delta
t=2\pi/\omega$, _i.e._ $\Delta=2\cos(2\pi A/\omega)$.

In the case when $q>\bar{r}$, we must identify the intervals $I_{i}$ at whose
end points $u(t;q)=\pm\bar{r}$. We start when $j$ is even. Owing to (17) we
have in the equation (29),

$\bar{f}_{j}(r)=\frac{1}{2}(mr+m\bar{r}),\quad\text{ if }\quad
r>\bar{r},\quad\bar{f}_{j}(r)=\frac{1}{2}(mr-m\bar{r}),\quad\text{ if }\quad
r<-\bar{r},$

therefore the two angular frequencies of the Hill equation (30) are

$A_{0}=\sqrt{\beta j^{2}+\gamma m},\quad\text{ if }\quad|u|\geq\bar{r},\qquad
A_{1}=\sqrt{\beta j^{2}+2\gamma m},\quad\text{ if }\quad|u|<\bar{r}.$

To determine the intervals in which the two cases occur, we note that
$\bar{f}_{j}(r)$ is an odd function, so that the potential for the equation
(29) is even. Therefore it is enough to compute $u(t;q)$ for a quarter of a
period, starting from $t=0$ to the first positive time $\,t_{1}=t_{1}(q)$ when
$\,u(t_{1};q)=0$, so that $\tau=4t_{1}$.

If we call $t_{0}=t_{0}(q)<t_{1}(q)$ the first positive time such that
$u(t_{0};q)=\bar{r}$, then the cycle $[-t_{0},\tau-t_{0}]$ is the union of
disjoint intervals $\cup_{i=0}^{3}I_{i}$, such that

$|u|\geq\bar{r},\quad t\in I_{0}\cup I_{2};\qquad|u|\leq\bar{r}\quad t\in
I_{1}\cup I_{3}.$

Their length are $\Delta t_{i}=2t_{0}$, for $i=0,2$, and $\Delta
t_{i}=2(t_{1}-t_{0})$, for $i=1,3$. As a consequence, the potential of the
equation (30) is given by

$Q(t)=A^{2}_{0}\mathbbm{1}_{I_{0}}(t)+A^{2}_{1}\mathbbm{1}_{I_{1}}(t)+A^{2}_{0}\mathbbm{1}_{I_{2}}(t)+A^{2}_{1}\mathbbm{1}_{I_{3}}(t),\qquad
t\in[-t_{0},\tau-t_{0}],$

and we have 4 transition matrices with $L_{0}=L_{2}$, $L_{1}=L_{3}$, so that

$M=(L_{1}\,L_{0})^{2},\qquad\Delta(q)={\rm tr}M.$

As for the computation of $t_{0}$, $t_{1}$, we define the 2 angular
frequencies of the equation (29),

$\omega_{0}=\sqrt{\alpha j^{4}+m},\qquad\omega_{1}=\sqrt{\alpha j^{4}+2m},$

so that $u$ satisfies the equation $\,\ddot{u}+\omega_{0}^{2}u=-m\bar{r}\,$ on
$I_{0}$, and $\ddot{u}+\omega_{1}^{2}u=0$ on $I_{1}$. On the first interval
$I_{0}=[-t_{0},t_{0}[$, we obtain

$u(t;q)=\left(q+\frac{m\bar{r}}{\omega_{0}^{2}}\right)\cos\omega_{0}t-\frac{m\bar{r}}{\omega_{0}^{2}}.$

By solving the equation $\,u(t;q)=\bar{r}$, we get

$t_{0}=\frac{1}{\omega_{0}}{\rm
acos}\,\left(\frac{\omega_{0}^{2}\bar{r}+m\bar{r}}{\omega_{0}^{2}q+m\bar{r}}\right).$

On the second interval $I_{1}$, $u$ solves the equation
$\ddot{u}+\omega_{1}^{2}u=0$, with initial data $u(t_{0};q)$ and
$\dot{u}(t_{0};q)$. Thus

$u(t;q)=\bar{r}\cos(\omega_{1}(t-t_{0}))+\dot{u}(t_{0};q)\sin(\omega_{1}(t-t_{0}))/\omega_{1}.$

Since we have

$\dot{u}(t_{0};q)=-\omega_{0}(q+\frac{m\bar{r}}{\omega_{0}^{2}})\sqrt{1-\cos^{2}\omega_{0}t_{0}}=-\sqrt{(q-\bar{r})[\omega_{0}^{2}(q+\bar{r})+2m\bar{r}]},$

by solving the equation $\dot{u}(t;q)=0$, we get

$t_{1}-t_{0}=\frac{1}{\omega_{1}}{\rm
atan}\,\left(\frac{\omega_{1}\bar{r}}{B(q)}\right),\qquad
B(q)=\sqrt{(q-\bar{r})[\omega_{0}^{2}(q+\bar{r})+2m\bar{r}]}.$

Now we come to the calculations in the case when $j$ is odd. The function
$\bar{f}_{j}(r)$ in (29), for $|r|>\bar{r}$, is given by

$\displaystyle\bar{f}_{j}(r)=\frac{m}{2}(1+1/j)r+\frac{m}{2}(1-i/j)\bar{r},\quad
r>\bar{r},$
$\displaystyle\bar{f}_{j}(r)=\frac{m}{2}(1-1/j)r-\frac{m}{2}(1+1/j)\bar{r},\quad
r<-\bar{r}.$

The procedure is similar to the previous case, only a bit longer, because
$\bar{f}_{j}(r)$ is no more odd, so that we need to compute $u(t;q)$ on a half
period. We define $t_{0}$ the first positive time such that
$u(t_{0};q)=\bar{r}$, and again we consider the cycle $[-t_{0},\tau-t_{0}]$
which is the union $\cup_{i=0}^{3}I_{i}$ of disjoint intervals, such that

$u\geq\bar{r},\quad t\in I_{0};\quad|u|\leq\bar{r},\quad t\in I_{1}\cup
I_{3};\quad u\leq-\bar{r},\quad t\in I_{2}.$

The Hill equation (30) is now a four steps Meissner equation with potential

$Q(t)=A^{2}_{0}\mathbbm{1}_{I_{0}}(t)+A^{2}_{1}\mathbbm{1}_{I_{1}}(t)+A^{2}_{2}\mathbbm{1}_{I_{2}}(t)+A^{2}_{1}\mathbbm{1}_{I_{3}}(t),\qquad
t\in[-t_{0},\tau-t_{0}],$

where

$A_{0}=\sqrt{\beta j^{2}+\gamma m(1+1/j)},\,\,A_{1}=A_{3}=\sqrt{\beta
j^{2}+2\gamma m},\,\,A_{2}=\sqrt{\beta j^{2}+\gamma m(1-1/j)},$

then the matrix $M$ becomes $M=L_{1}\,L_{2}\,L_{1}\,L_{0}$.

We define the 3 angular frequencies for the equation (29),

$\omega_{0}=\sqrt{\alpha j^{4}+m(1+1/j)},\quad\omega_{1}=\sqrt{\alpha
j^{4}+2m},\quad\omega_{2}=\sqrt{\alpha j^{4}+m(1-1/j)},$

and the constants

$D=m(1-1/j)\bar{r},\quad E=m(1+1/j)\bar{r},$

so that $u$ satisfies the equation $\ddot{u}+\omega_{0}^{2}u=-D$ on $I_{0}$,
the equation $\ddot{u}+\omega_{1}^{2}u=0$ on $I_{1}\cup I_{3}$, and
$\ddot{u}+\omega_{2}^{2}u=E$ on $I_{2}$.

Proceeding as above, for $t\in I_{0}$ we have,

$u(t;q)=\left(q+\frac{D}{\omega_{0}^{2}}\right)\cos\omega_{0}t-\frac{D}{\omega_{0}^{2}},$

thus, by finding $t_{0}$, we obtain half the length of $I_{0}$, that is

$\frac{\Delta t_{0}}{2}=\frac{1}{\omega_{0}}{\rm
acos}\,\left(\frac{\omega_{0}^{2}\bar{r}+D}{\omega_{0}^{2}q+D}\right).$

For $t\in I_{1}$, starting with the initial conditions at time $t_{0}$, we get
the expression

$u(t;q)=\bar{r}\cos(\omega_{1}(t-t_{0}))-\frac{B}{\omega_{1}}\sin(\omega_{1}(t-t_{0})),\quad
B=-u(t_{0};q)=\sqrt{(q-\bar{r})[\omega_{0}^{2}(q+\bar{r})+2D]}.$

Thus, by finding the first time greater than $t_{0}$ such that $u(t;q)=0$, we
obtain half the length of $I_{1}$:

$\frac{\Delta t_{1}}{2}=\frac{1}{\omega_{1}}{\rm
atan}\,\left(\frac{\omega_{1}\bar{r}}{B}\right).$

Finally, for $t\in I_{2}$, we need to solve the third equation
$\ddot{u}+\omega_{2}^{2}u=E$ with initial data at $t^{*}=t_{0}+\Delta t_{1}$,
that is $u(t^{*};q)=-\bar{r}$, $\dot{u}(t^{*};q)=\dot{u}(t_{0};q)=-B$.
Therefore, we have

$u(t;q)=-\left(\bar{r}+\frac{E}{\omega_{2}^{2}}\right)\cos\omega_{2}(t-t^{*})-\frac{B}{\omega_{2}}\sin\omega_{2}(t-t^{*})+\frac{E}{\omega_{2}^{2}}.$

By finding the first time greater than $t^{*}$ such that $\dot{u}(t;q)=0$, we
obtain half the length of $I_{2}$:

$\frac{\Delta t_{2}}{2}=\frac{1}{\omega_{2}}{\rm
atan}\,\left(\frac{\omega_{2}B}{\omega_{2}^{2}\bar{r}+E}\right).$

## References

  * [1] Ammann O.H., von Kármán T., Woodroff G.B., The Failure of the Tacoma Narrow Bridge, Tecnical Report,Federal Works Agency, Washington, DC, 1941, https://resolver.caltech.edu/CaltechAUTHORS:20140512-105559175
  * [2] Arioli G., Gazzola F., A new mathematical explanation of what triggered the catastrophic torsional mode of the Tacoma Narrows Bridge, Appl. Math. Model. 39 (2015), no. 2, 901-912, https://doi.org/10.1016/j.apm.2014.06.022
  * [3] Arnol’d, V. I. Mathematical methods of classical mechanics. Graduate Texts in Mathematics, 60. Springer-Verlag, New York-Heidelberg, 1978, https://doi.org/10.1007/978-1-4757-1693-1
  * [4] Berchio E., Gazzola F., A qualitative explanation of the origin of torsional instability in suspension bridges, Nonlinear Anal. 121 (2015), 54-72, https://doi.org/10.1016/j.na.2014.10.026
  * [5] Brown B.M., Eastham M.S.P., Schmidt K.M., Periodic differential operators. Operator Theory: Advances and Applications, 230. Birkhäuser/Springer, Basel, 2013, https://doi.org/10.1007/978-3-0348-0528-5
  * [6] Cazenave T., Weissler F.B., Unstable simple modes of the nonlinear string, Quart. Appl. Math. 54 (1996), no. 2. 287-305, https://doi.org/10.1090/qam/1388017
  * [7] Eastham M.S.P., The spectral theory of periodic differential equations. Texts in Mathematics. Scottish Academic Press, Edinburgh; Hafner Press, New York, 1973
  * [8] Falocchi A., Torsional instability in a nonlinear isolated model for suspension bridges with fixed cables and extensible hangers , IMA J. Appl. Math. 83 (2018), no. 6, 1007-1036, https://doi.org/10.1093/imamat/hxy032
  * [9] Falocchi A., Torsional instability and sensitivity analysis in a suspension bridge model related to the Melan equation, Commun. Nonlinear Sci. Numer. Simul. 67 (2019), 60-75, https://doi.org/10.1016/j.cnsns.2018.07.005
  * [10] Gazzola F., Mathematical Models for Suspension Bridges: Nonlinear Structural Instability, MS&A. Modeling, Simulation and Applications, 15. Springer, Cham, 2015, https://doi.org/10.1007/978-3-319-15434-3
  * [11] Hirsch M.W., Smale S., Differential equations, dynamical systems, and linear algebra. Pure and Applied Mathematics, Vol. 60. Academic Press, New York-London, 1974
  * [12] Hochstadt H., A Special Hill’s Equation with Discontinuous Coefficients, N.Y.U. Inst. Math. Sci., Div. of EM Res., Report BR-32 (1960)
  * [13] Loud W., Stability regions for Hill’s equation, J. Differential Equations, 19 (1975), 226-241, https://doi.org/10.1016/0022-0396(75)90003-0
  * [14] Magnus W., Winkler S., Hill’s equation. Interscience Tracts in Pure and Applied Mathematics, No. 20 Interscience Publishers John Wiley & Sons, New York-London-Sydney 1966,
  * [15] Marchionna C., Panizzi S., An instability result in the theory of suspension bridges, Nonlinear Anal. 140 (2016), 12-28, https://doi.org/10.1016/j.na.2016.03.003
  * [16] Marchionna C., Panizzi S., An instability result in the theory of suspension bridges, Integral methods in science and engineering. Vol. 1. Theoretical techniques, 193-203, Birkhäuser/Springer, Cham, 2017, https://doi.org/10.1007/978-3-319-59384-5_17
  * [17] Marchionna C., Panizzi S., On the instability tongues of the Hill equation coupled with a conservative nonlinear oscillator, J. Math. Anal. Appl. 479, (2019), no. 2, 2139 - 2164, https://doi.org/10.1016/j.jmaa.2019.07.045
  * [18] McKenna P.J., Walter W., Nonlinear oscillations in a suspension bridge, Arch. Ration. Mech. Anal., 98 (1987), 167-177
  * [19] McKenna P.J., Tuama C. Ó, Large torsional oscillations in suspension bridges visited again: vertical forcing creates torsional response, Amer. Math. Monthly 108 (2001. https://doi.org/10.1080/00029890.2001.11919805
  * [20] McLachlan, N. W., Theory and application of Mathieu functions. Dover Publications, Inc., New York 1964
  * [21] Moore K.S., Large torsional oscillations in a suspension bridge: multiple periodic solutions to a nonlinear wave equation, SIAM J. Math. Anal. 33 (2002), no. 6, 1411-1429, https://doi.org/10.1137/S0036141001388099
  * [22] Opial Z., Continuous parameter dependence in linear systems of differential equations, J. Differential Equations 3 (1967), 571-579, https://doi.org/10.1016/0022-0396(67)90017-4
  * [23] Šermr J., _On differentiation of a Lebesgue integral with respect to a parameter_ , Math. Appl. (Brno) 1 (2012), no. 1, 91-116, https://doi.org/10.13164/ma.2012.06
  * [24] Gan Shaobo, Zhang Meirong, Resonance Pockets of Hill’s Equations with Two-Step Potentials, SIAM J. Math. Anal., 32(3), (2000), 651–664, https://doi.org/10.1137/S0036141099356842
  * [25] Van Der Pol B., Strutt M. J. O., On the stability of the solutions of Mathieu’s equation, The London, Edinburgh and Dublin Phil. Mag. 7th Series 5 (1928), 18-38, https://doi.org/10.1080/14786440108564441
  * [26] Weinstein M. I. and Keller J. B., Hill’s equation with a large potential, SIAM J. Appl. Math. 45 (1985), no.2, 200-214, https://doi.org/10.1137/0145011