MATLAB TUTORIAL, Part 1.2: Linear Equations

Prof. Vladimir A. Dobrushkin

This tutorial contains many matlab scripts.
You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to

Email Vladimir Dobrushkin

Linear Equations

A linear first order ordinary differential equation (ODE) can be used as a mathematical model for a variety of phenomena, either physical or non-physical. Examples of such phenomena include the following: heat flow problems (thermodynamics), simple electrical circuits (electrical engineering), force problems (mechanics), rate of bacterial growth (biological science), rate of decomposition of radioactive material (atomic physics), crystallization rate of a chemical compound (chemistry), and rate of population growth (statistics). Most of these problems, however, appear in systems of differential equations that are considered in the second course.

A differential equation, written in the normal form:

\[ y' + a(x) y = f(x) , \]

is called the linear equation, where a(x) is a coefficient function and f(x) is the forcing, driving, input, or nonhomogeneous term. If the driving term f(x) is not identically zero, then the linear equation is called nonhomogeneous/inhomogeneous or driven. Otherwise, it is called the homogeneous equation.

Theorem: Let a(x) and f(x) be continuous functions on the open interval (a,b), and let \( x_0 \in (a,b) . \) Then for each \( x \in (a,b) \) there exists a unique solution \( y = \phi (x) \) to the differential equation \( y' + a(x)\,y = f(x) \) that also satisfies the initial value condition that \( y(x_0 ) = y_0 \) for any real number y0. Moreover, this initial value problem has no singular solution. ■

Note that if the interval in the above theorem is the largest possible interval on which a(x) and f(x) are continuous, then the interval is the interval of validity for the solution. This means, that for linear first order differential equations, we won't need to actually solve the differential equation in order to find the largest possible interval where the solution exists and continuous (such interval is called the validity interval). Notice as well that the interval of validity will depend only partially on the initial condition. The interval must contain x0, but the value of y0 has no effect on the interval of validity.

Example: validity interval: Determine the validity interval for the initial value problem

\[ \left( x^2 -4 \right) y' + 3\,y = \ln |7-x| , \qquad y(4) = 5. \]

Solution: First, in order to use the theorem to find the interval of validity, we must rewrite the differential equation in the normal form given in the theorem. So we need to divide out by the coefficient of the derivative.

\[ y' + \frac{3}{\left( x^2 -4 \right)}\,y = \frac{\ln |7-x|}{\left( x^2 -4 \right)} , \qquad y(4) = 5. \]

Next, we need to identify where the two functions \( a(x) = \frac{3}{\left( x^2 -4 \right)} \) and \( f(x) = \frac{\ln |7-x|}{\left( x^2 -4 \right)} \) are not continuous. This will allow us to find all possible intervals of validity for the differential equation. So, a(x) is discontinuous at \( x= \pm 2 , \) while f(x) is also undefined at x = 7. Now, with these points in hand, we can break up the real number line into four intervals where both a(x) and f(x) will be continuous. These four intervals are,

\[ -\infty < x < -2, \quad -2 < x < 2, \quad 2 < x < 7, \quad 7 < x < \infty . \]
The endpoints of each of the intervals are points where at least one of the two functions is discontinuous. This will guarantee that both functions are continuous everywhere in each interval. Finally, let's identify the actual interval of validity for the given initial value problem. The actual interval of validity is the interval that will contain x0 = 4. So, the interval of validity for the initial value problem is
\[ 2 < x < 7 . \]
In this last example we need to be careful to not jump to the conclusion that the other three intervals cannot be intervals of validity. By changing the initial condition, in particular the value of x0, we can make any of the four intervals the interval of validity. ■

A linear homogeneous (also called not driven) equation
\[ y' + a(x) y = 0 \]
is also a separable equation; its solution can be obtained by separation of variables:
\[ \frac{{\text d}y}{y} = -a(x)\,{\text d}x \qquad \Longrightarrow \qquad y(x) = C\,\exp \left\{ -\int a(x)\,{\text d}x \right\} . \]
where C is an arbitrary constant.

Example: Consider the homogeneous linear equation

\[ y' (x) -5\, y =0 . \]
Separating variables, we obtain
\[ \frac{{\text d}y}{y} = 5\,{\text d}x \qquad \Longrightarrow \qquad y(x) = C\,e^{5x} . \]
We implement its solution using the following Mathematica commands.

ode[x_,y_] = (y'[x] - 5 y[x] == 0)
SolRule = DSolve[ode[x, y], y[x], x]
yh[x_] = Simplify[y[x] /. SolRule[[1]]]            (* solution of the homogeneous eq *)
Simplify[ode[x, yh]]                                       (* to check the answer *)

For nonhomogeneous linear equation, there are known two systematic methods to find their solutions: integrating factor method and the Bernoulli method.

Integrating factor method allows us to reduce a linear differential equation in normal form \( y' + a(x)\,y = f(x) \) to an exact equation. There always exists an integrating factor μ(x) as a function of x:

\[ \mu (x) = \exp \left\{ \int a(x)\,{\text d}x \right\} , \]
which is a solution of the first order homogeneous linear equation:
\[ \frac{{\text d}\mu (x)}{{\text d}x} - a(x)\, \mu (x) =0 \qquad \Longrightarrow \qquad \frac{{\text d}\mu (x)}{\mu} = a(x)\, {\text d}x. \]
Multiplying both sides of the nonhomogeneous differential equation \( y' + a(x)\,y = f(x) \) by the integrating factor μ(x) yields
\[ \mu \,\frac{{\text d}y (x)}{{\text d}x} + a(x)\, \mu (x)\, y = \mu (x)\, f(x) \qquad \Longleftrightarrow \qquad \frac{{\text d}}{{\text d}x} \left( \mu (x)\, y(x) \right) = \mu (x)\, f(x) , \]
which is an exact equation. Integration yields the general solution
\[ \mu (x)\, y(x) = \int f(x)\, \mu (x) \,{\text d}x +C \qquad \Longrightarrow \qquad y(x) = \frac{1}{\mu (x)} \, \int f(x)(x)\, \mu (x) \,{\text d}x + \frac{C}{\mu (x)} . \]

The Bernoulli method suggests to seek a solution of the inhomogeneous linear differential equation (it does not matter whether the equation is in normal form or not) \( y' + a(x)\,y = f(x) \) in the form of the product of two functions:

\[ y(x) = u(x)\, v(x) , \]
where u(x) is a solution of its homogeneous part:
\[ u' + a(x)\, u =0 , \]
and v(x) is the general solution of the separable and exact equation:
\[ u\, v' = f(x) \qquad \Longrightarrow \qquad \frac{{\text d}v}{{\text d}x} = u^{-1} f(x) . \]
Integrating the latter and multiplying it by u, we obtain the general solution:
\[ y(x) = u(x)\, \int u^{-1} f(x) \,{\text d}x + C\, u(x) , \]
where C is a constant of integration.

Example: Consider the nonhomogeneous equation

\[ y' (x) -5\, y =27\,x^3 e^{2x} . \]
An integrating factor μ(x) for this differential equation is a solution of the separable equation
\[ \mu' (x) +5\, \mu (x) =0 \qquad \Longrightarrow \qquad \frac{{\text d}\mu}{\mu} = - 5\,{\text d}x . \]
Since we need only one integrating factor to solve the given differential equation, we can choose it as \( \mu (x) = e^{-5x} . \) Upon multiplying by it both sides, we reduce the given equation to an exact equation:
\[ e^{-5x} \,\frac{{\text d}y (x)}{{\text d}x} -5\, e^{-5x}\, y = e^{5x}\, 27\,x^3 e^{2x} \qquad \Longleftrightarrow \qquad \frac{{\text d}}{{\text d}x} \left( e^{-5x}\, y(x) \right) = 27\,x^3 e^{-3x} , \]
Integrating, we obtain the general solution, which we also get by executing the Mathematica codes:
nde[x_, y_] = (y'[x] - 5 y[x] == 27 x^3 Exp[2 x])
SolnRule = DSolve[nde[x, y], y[x], x]
yn[x_] = Simplify[y[x] /. SolnRule[[1]]]
Out[11]= -E^(2 x) (2 + 6 x + 9 x^2 + 9 x^3) + E^(5 x) C[1]

Its solution can be obtained in one line Mathematica code:

yn[x_] = Simplify[(y[x] /. DSolve[nde[x, y], y[x], x][[1]])]

Now we apply the Bernoulli method: \( y(x) = u(x)\,v(x) , \) where u is a solution of the homogeneous equation

\[ u' (x) -5\, u = 0 \qquad \Longrightarrow \qquad u(x) = e^{5x} . \]
Then for v(x) we get an exact equation:
\[ u\,v' (x) = 27\,x^3 e^{2x} \qquad \Longrightarrow \qquad v'(x) = 27\,x^3 e^{-3x} . \]
Simple integration
Integrate[27*x^3 * Exp[-3*x],x]
yields the function v(x), upon multiplying it by u(x), we get the general solution:
\[ y(x) = C\, e^{5x} + e^{2x} \left( 2+ 6x + 9\, x^2 + 9\, x^3 \right) . \]

Example: Consider another example for a linear equation with variable coefficients:

\[ x\, y' -5\,y = 27\,x^7 e^{x} . \]
Of course, we can find its solution using a standard command DSolve:
SolxRule = DSolve[x y'[x] - 5 y[x] == 27 x^7 * Exp[x], y[x], x]
y3[x_] = Simplify[y[x]/.SolxRule[[1]]]
Out[11]= x^5 (27 E^x (-1 + x) + C[1])

Checking the answer:

Simplify[x y3'[x] -5 y3[x] == 27 x^7 * Exp[x]]
we conclude that the obtained function is indeed its general solution.

We use an integrating factor method by solving the corresponding homogeneous/separable equation for μ(x):

\[ x\, \mu' +5\,\mu = 0 \qquad \Longrightarrow \qquad \frac{{\text d}\mu}{\mu} = - 5\frac{{\text d}x}{x} \qquad \Longrightarrow \qquad \mu (x) = x^{-5} . \]
To obtain the exact equation, we need first to reduce the given equation to the normal form by isolating the derivative: \( y' - 5\,x^{-1}\,y = 27\,x^6 e^x . \) Multiplying both sides of the latter by μ(x), we obtain an exact equation:
\[ \frac{\text d}{{\text d}x} \left[ \mu (x)\,y(x) \right] = \mu (x)\,27\,x^6 e^{x} \qquad \Longleftrightarrow \qquad \frac{\text d}{{\text d}x} \left[ x^{-5}\,y(x) \right] = 27\,x\, e^{x} . \]
Integration yields
Integrate[27*x * Exp[x],x]
27 E^x (- 1 + x )
Multiplying the above function by x5, we obtain the general solution:
\[ y(x) = 27\, x^{5} e^x \left( -1 + x \right) + C\,x^{5} , \]
where C is an arbitrary constant.

Now we demonstrate the Bernoulli method: \( y(x) = u(x)\, v(x) , \) where u(x) is a solution of the homogeneous part \( x\, u' -5\,u =0 \) and v(x) is obtained upon solving the separable equation \( x\,u\, v' = 27\, x^7 \, e^x . \) Integrating these sequential equations

\[ \begin{split} x\, u' &= 5\, u \qquad \Longrightarrow \qquad u = x^5 , \\ x\,u\, v' &= 27\, x^7 \, e^x \qquad \Longleftrightarrow \qquad v' = 27\, x\, e^x \qquad \Longrightarrow \qquad v= 27\, e^x \left( x-1 \right) , \end{split} \]
we find functions u and v; their product gives us the required solution:
\[ y(x) = 27\, x^{5} e^x \left( -1 + x \right) + C\,x^{5} . \]

Example: Solve the IVP: \( 3ty' +2y=t^2 , \quad y(1)=1 . \) We try to find its solution using Mathematica:

Clear[t, y];
t0 = 1; y0 = 1;
p[t_] = 2/3/t; q[t_] = t/3;
P[t_] = Exp[Integrate[p[s], {s, t0, t}]];
y[t_] = 1/P[t]*(t0 + Integrate[P[u]*q[u], {u, t0, t}])
Out[2]= ConditionalExpression[ConditionalExpression[(1 + 1/8 (-1 + t^(8/3)))/t^(2/3), Re[t] >= 0 || t \[NotElement] Reals],
(Re[t] >= 0 || t \[NotElement] Reals) && (Re[u] >= 0 || u \[NotElement] Reals)]

Another option is to apply DSolve directly:
DSolve[{3*t*y'[t]+2*y[t] == t^2 , y[1]==1},y,t]
which gives
\[ y(t) = \frac{7 + t^{8/3}}{8\,t^{2/3}} + C\, t^{-2/3} . \]

Example: Using Mathematica, solve the linear differential equation: \( y' = e^{-2x} -3\,y , \) and plot some its solutions.

Clear[x,y,gensol]                                                (* clear all variables *)
gensol=DSolve[y'[x]==Exp[-2 x]-3 y[x],y[x],x]
Out[11]= {{y[x] -> E^(-2 x) + E^(-3 x) C[1]}}
gensol[[1,1,2]]
Out[12]= E^(-2 x) + E^(-3 x) C[1]

In DSolve command, the first argument (y'[x]==Exp[-2 x]-3 y[x]) represents the differential equation, the second argument (y[x]) instructs Mathematica that we are solving for      y=y(x),      and the third argument (x) instructs Mathematica that the independent variable is x.
Note that gensol is a nested list. The first part of gensol, extracted with gensol[[1]], is the list (y(x)->E^(-2 x) + E^(-3 x) C[1]);
and the first part of this list, extracted with gensol[[1,1,1]], is y(x) while the second part of this list
(which represents the formula for the solution), extracted with gensol[[1,1,2]], is y=E^(-2 x) + E^(-3 x) C[1]

toplot=Table[(gensol[[1,1,2]])/.C[1]->i,{i,-2,2,0.25}];
pp=Plot[Evaluate[toplot],{x,-1/2,1},PlotRange->{-3/4,3/4},AspectRatio->1,
AxesStyle->GrayLevel[0.5],PlotStyle->{{GrayLevel[0.4],Thickness[0.01]}}, AxesLabel->{x,y}]

pd = VectorPlot[{1, Exp[-2 x] - 3 y}, {x, -1/2, 1}, {y, -3/4, 3/4},
Axes -> Automatic, VectorScale -> (1 &), AxesOrigin -> {0, 0}, VectorStyle -> Arrowheads[0]]

Show[pd, pp]


Example:In mining, “mine tailing” are what is left after everything valuable (such as a mineral, coal, or oil) has been removed. The material that is left over after the minerals, coal, or oil is extracted often presents a great deal of hazard to the environment. One of the ways of processing mine tailing is to store them in a pong; this method is commonly used when water is used in the mining extraction. This method allows any particles that are suspended in the water to settle at the bottom of the pond. The water can then be treated and recycled.

Suppose that we have a gold mining operation and we are storing our tailing in a pond that has an initial volume of 100,000 cubic meters. When we begin our operation, the pond is filled with clean water. The pond has a stream flowing into it, and water is also being pumped out of the pond. Chemicals are used as a way to process gold ore, which is the material being extracted in this operation. Chemicals that are used, like sodium cyanide, are often highly poisonous and harmful to the environment. Thus, the water must be treated before it is released into the watershed. Suppose that 3,000 cubic meters per day flow into the pond from stream, and 3,000 cubic meters are pumped out of the pond each day to be processed and recycled. Since inflow and outflow are equal, the water level of the pond remains constant.

At time t=0, the water from stream becomes contaminated with chemicals from the mining operation at a rate of 10 kilograms of chemicals per 1000 cubic meters. We will assume that water in our tailing pond is well mixed so that the concentration of chemicals throughout the pond is uniform. In addition, any matter pumped into the pond from the stream settles to the bottom of the pond at a rate of 100 cubic meters per day. Thus, the volume of the pond is reduced by 100 cubic meters each day, and will become full after 500 days of operation. We shall assume that the particulate matter and the chemicals are included in the 1000 cubic meters that flow into the pond from the stream each day.

We want to find a differential equation that will model the amount of chemicals in the tailing pond at any particular time. Let y(t) be the amount of chemicals in the pond at time t. Then \( {\text d}y/{\text d}t \) is the difference between the rate at which the chemicals enter the pond and the rate at which the chemicals leave the pond.

\[ \frac{{\text d} y}{{\text d} t} = \mbox{rate in} - \mbox{rate out} . \]

Since water flows into the pond from the stream at a rate of 1000 cubic meters per day, the rate as which the chemicals enter the pond is 10 kilograms per day. On the other hand, the rate at which the chemicals leave the pond will depend on the amount of chemicals in the pond at time t. The volume of the pond is decreasing due to sediment, and at time t it is V(t)=100000−100t. Thus, the concentration of chemicals in the pond at time t is y/(100000−100t), and the rate at which the chemicals are flowing out of the pond to be recycled is

\[ 1000 \left( \frac{y}{100000-100t} \right) = \frac{10\,y}{1000 -t} . \]
Hence, the differential equation that models the amount of chemical in the tailing pond at time t is
\[ \frac{{\text d} y}{{\text d} t} = 10-\frac{10\,y}{1000 -t} . \]
Of course, we will have to cease mining operations once the pond is full, since there will only be water in the pond if \( V(t)=100000−100t \ge 0; \) that is, when \( 0 \le t< < 1000. \)

Notice that the above equation is not autonomous. In fact, it is not even separable. We will have to use a different approach to find a solution. First, we will rewrite the equation in the form suitable for an integrating factor

\[ \frac{{\text d} y}{{\text d} t} + \frac{10\,y}{1000 -t} = 10 . \]
If we multiply both sides of this equation by (1000−t)-10, we obtain
\[ \left( 1000 -t \right)^{-10} \frac{{\text d} y}{{\text d} t} + 10\,y \left( 1000 -t \right)^{-11} = 10 \left( 1000 -t \right)^{-10} . \]
We now make the crucial observation that the product rule applies to the left-hand side of our equation,
\[ \left( 1000 -t \right)^{-10} \frac{{\text d} y}{{\text d} t} + 10\,y \left( 1000 -t \right)^{-11} = \frac{{\text d}}{{\text d}t} \left[ \left( 1000 -t \right)^{-10} y(t) \right] . \]
Thus, our equation becomes exact and integrating both sides, we obtain
\[ \left( 1000 -t \right)^{-10} y(t) = 10 \int {\text d}t \left( 1000 -t \right)^{-10} = \frac{10}{9} \left( 1000 -t \right)^{-9} +C, \]
where C is an arbitrary constant. Solving for y, we obtain
\[ y(t) = \frac{10}{9} \left( 1000 -t \right) +C \left( 1000 -t \right)^{10} . \]
Since y(0)=0, we can quickly determine that \( C=−(10/9)\, 1000^{−9} \) and that the solution to our initial value problem is
\[ y(t) = \frac{10}{9} \left( 1000 -t \right) \left[ 1 - \left( \frac{1000-t}{1000} \right)^{9} \right] . \]
The graph of the solution to our differential equation (Figure ) fits the situation. Initially, there are no chemicals in the pond, but y(t) quickly increases. However, the amount of chemicals decreases as the pond begins to fill with sediment. Eventually, there are no chemicals at t=1000.
f[t_] = 10/9*(1000 - t)*(1 - (1 - t/1000)^9)
Plot[f[t], {t, 0, 1009}, PlotStyle -> Thick]
Mine tailing.

Example: Mixing Models. Many applications involve the mixing of two or more substances together. We can model how petroleum products are mixed together in a refinery, how various ingredients are mixed together in a brewery, or how greenhouse gases mix and move across various layers of the earth's atmosphere. Basically, it consists of finding a formula for the amount of some "pollutant" in a container, into which pollutant is entering at a fixed rate and also flowing out at a fixed rate. The general physical rule used to describe this situation is

\[ \mbox{rate of change of pollutant per unit time} = \mbox{rate in} - \mbox{rate out} . \]
If we denote by y(t) the amount of pollutant in the container at time t, its rate of change per unit time is given by \( \frac{{\text d}y}{{\text d}t} = \dot{y} ; \) therefore, the problem leads to a differential equation for the function y(t).

Suppose that a 400-liter tank initially contains 200 liters of salt water containing 2 kilograms of salt. A brine mixture containing 1/10 kilograms of salt per liter flows into the top of the tank at a rate of 4 liters per minute. A well-mixed solution leaves the tank at rate of 3 liters per minute. We wish to know how much salt is in the tank, when the tank is full.

To construct our model, we will let t be the time (measured in minutes) and set up a differential equation that will measure how fast the amount of salt at time t, y(t), is changing. We have the initial condition y(0) = 5, and

\begin{eqnarray*} \frac{{\text d} y}{{\text d} t} &=& \mbox{amount of salt in} - \mbox{amount of salt out} \\ &=& \frac{4}{10} - 3\,\frac{y}{V(t)} \end{eqnarray*} where V(t) is the volume at time t. The expression x/V(t) is the amount of salt in one liter at time t. We have V(t)=200+t, since the tank starts with 200 liters and four liters are pumped into the tank per minute while three liters leaves the tank during the same time interval. Thus, our differential equation becomes
\[ \frac{{\text d} y}{{\text d} t} = \frac{4}{10} - \frac{3}{200+t}\,y . \]
Using the Bernoulli method, we seek its solution as y = uv, where u is a solution of homogeneous equation
\[ \frac{{\text d} u}{{\text d} t} = - \frac{3}{200+t}\,u \qquad \Longrightarrow \qquad \frac{{\text d} u}{u} = - \frac{3}{200+t}\, {\text d} t . \]
Integration yields
\[ \ln u = - 3\,\ln \left( 200+t \right) = \ln \left( 200+t \right)^{-3} \qquad \Longrightarrow \qquad u(t) = \left( 200+t \right)^{-3} . \]
Then for v(t), we get the separable equation
\[ u\, \frac{{\text d} v}{{\text d} t} = \frac{4}{10} \qquad \Longrightarrow \qquad \frac{{\text d} v}{{\text d} t} = \frac{4}{10} \left( 200+t \right)^{3} \qquad \Longrightarrow \qquad v(t) = \frac{1}{10} \left( 200+t \right)^{4} +C , \]
where C is an arbitrary constant. Therefore, the general solution becomes
\[ y(t) = u(t)\,v(t) = \frac{1}{10} \left( 200+t \right) + C \left( 200+t \right)^{-3} . \]
Our initial condition, x(0)=2 tells us that C=-18*2003 and
\[ y(t) = u(t)\,v(t) = \frac{1}{10} \left( 200+t \right) - \frac{144000000}{\left( 200+t \right)^{3}} . \]
Thus, when the tank is full, t=200 and the amount of salt in the tank is \( y(200)= \frac{151}{4} = 37.75 \) kilograms.

Example: gas in magma affects volcanic eruptions. Many andesitic volcanoes exhibit effusive eruption activity, with magma volumes as large as 107 -- 109 m3 erupted at rates of 1 -- 10 m3/s over periods of years or decades. During such eruptions, many complex cycles in eruption rates have been observed, with periods ranging from hours to years. Longerterm trends have also been observed, and are thought to be associated with the continuing recharge of magma from deep in the crust and with waning of overpressure in the magma reservoir. Here we present a model which incorporates effects due to compressibility of gas in magma. The eruption duration and volume of erupted magma may increase by up to two orders of magnitude if the stored internal energy associated with dissolved volatiles can be released into the magma chamber. This mechanism would be favored in shallow chambers or volatile-rich magmas and the cooling of magma by country rock may enhance this release of energy, leading to substantial increases in eruption rate and duration.

Consider a magma with bulk density ρ in a chamber of volume V undergoing a mass recharge rate Qi and eruption rate Q0. Conservation of mass indicates that

\[ \frac{{\text d}}{{\text d}t} \left( \rho\, V \right) = \rho\,\frac{{\text d} V}{{\text d}t} + V\, \frac{{\text d}\rho}{{\text d}t} = Q , \]
where \( Q = Q_i - Q_0 . \) The density of the magma, which consists of melt, crystals, and gas, can be written, in general form, as \( \rho = \rho \left[ p, T, x(T), N \right] , \) where p is pressure, T is temperature, x(T) is the mass fraction of crystals in the magma and N is the total mass fraction of volatiles. Depending on the pressure, a fraction of these volatiles are dissolved in the magma and the remainder exist in the gaseous phase. Differentiating this expression for the density and incorporating the result into the conservation of energy equation, we obtain:
\[ \frac{{\text d}V}{{\text d}t} + \frac{V}{\rho} \, \frac{\partial \rho}{\partial p} \, \frac{{\text d}p}{{\text d}t} = \frac{Q}{\rho} - \frac{V}{\rho} \,\frac{\partial \rho}{\partial T} \,\frac{{\text d} T}{{\text d}t} , \]
where the second term on the right-hand side represents the rate of change of volume associated with the change in density with temperature of the magma, which includes the effects of the change of crystal content with temperature, leading to an effective thermal expansion.

It is known that the rate of change in volume of the chamber, \( {\text d}V / {\text d}t , \) is related to the associated rate of change in pressure, \( {\text d}p / {\text d}t , \) as a result of deformation of the surrounding rock by:

\[ \frac{1}{\beta_r} \, \frac{{\text d}\rho}{{\text d}t} = \frac{1}{V} \, \frac{{\text d} V}{{\text d}t} , \]
where βr is the effective bulk modulus of the surrounding wall rock. The exact value of βr depends on the density of microfractures in the rock, but a value of 1010 Pa is typical. The equation describing this is
\[ \frac{{\text d} p}{{\text d}t} = \frac{\overline{\beta}}{V} \left( \frac{Q_i - \lambda p}{\rho} - \frac{V}{\rho} \, \frac{\partial \rho}{\partial T} \, \frac{{\text d} T}{{\text d}t} \right) , \]
where p is the chamber overpressure which is proportional to the volume eruption rate, Qi is the constant mass recharge rate, \( \overline{\beta} \) is the bulk modulus of the magma, V is the chamber volume, ρ is the density, and λ is is just a constant based on magma viscosity and the properties of the magma chamber, both of which do not change with time.

Although the initial equation has partial derivatives, it can be assumed for the sake of this model that temperature is constant throughout the eruption thus the equation can be simplified to an ODE and solved both numerically and analytically.

\[ \frac{{\text d} p}{{\text d}t} = \frac{\overline{\beta}}{V} \left( \frac{Q_i - \lambda p}{\rho} \right) . \]
Using separation of variables, we find its solution:
\[ p(t) = \left( p_0 - \frac{Q_i}{\lambda} \right) e^{- \lambda \overline{\beta} \,t/(\rho V)} + \frac{Q_i}{\lambda} . \]
This example is based on the article:
Herbert E. Huppert & Andrew W. Woods.,
The role of volatiles in magma chamber dynamics,
nature, Vol 420, 493--495, 2002.

%Constants throughout problem
x=0.4; %crystal content by mass percent
R= 462; %J/(K*Kg)
N=.02; %volatile constant by mass percent
V= 10; %km^3 chamber volume
rho_c= 2600; %kg/m^3 crystal density
rho_m= 2300; %kg/m^3 melt density
rho_r= 2500; %kg/m^3 rock above chamber density
H= 5000; % m depth of chamber
T= 1073; % K temperature of system
Beta_r= 10e10; % Pa bulk modulus of surrounding rock
s= 4e-6; %Pa^-.5
mu= 10e7; %Pa dynamic viscosity
S=0.1; %shape factor
r= 15;%m conduit radius What is a reasonable value for this???
A= pi*r^2;%m^2 conduit cross section
g= 9.81;%m/s^2 
p_b= rho_r * g* H; %background pressure of chamber
Qi=0; %mass recharge rate
N_iters=1000;
dt=.01;
%pre-allocating answer vectors
delta_p= zeros(N_iters+1,1); %over pressure
p=zeros(N_iters+1,1); %total pressure
rho_g= zeros(N_iters+1,1); %gas density based on ideal gas law
n=zeros(N_iters+1,1);%exsolved volatiles
rho=zeros(N_iters+1,1);%total density
lambda=zeros(N_iters+1,1);%constant relating over pressure to mass eruption rate
Beta=zeros(N_iters+1,1);%bulk modulus of compressible magma
Qo=zeros(N_iters+1,1);%mass outflow
%initial conditions
delta_p(1)= 10e6; %Pa initial condition
for i=1:N_iters
%values that depend on time
    p(i) = delta_p(i) + p_b;
    rho_g(i)=p(i)/(R*T);
    n(i)= N-s*p(i)^.5 *(1-x);
    if s*p(i)^.5 *(1-x)>=N 
        n(i)=0;
    end 
    rho(i)=((n(i)/rho_g(i))+(1-n(i))*((x/rho_c)+((1-x)/rho_m)))^-1;    
    Beta(i)=((1/Beta_r)+((N*R*T*p(i)^(-2)-.5*R*s*T*(1-x)*p(i)^(-1.5)+.5*(x/rho_c + (1-x)/ rho_m)*s*(1-x)*p(i)^(-.5))/(N*R*T*p(i)^(-1)-R*s*T*(1-x)*p(i)^(-.5)+(x/rho_c + (1-x)/ rho_m)-(x/rho_c + (1-x)/ rho_m)*N+(x/rho_c + (1-x)/ rho_m)*s*(1-x)*p(i)^.5)^2)/rho(i))^-1;
    lambda(i)=(rho(i)*S*A^2)/(H*mu);
    K1=(Qi-lambda(i)*delta_p(i))*(Beta(i)/(rho(i)*V))*dt;
        delta_p_rk=delta_p(i)+K1;
        p_rk = (delta_p(i)+K1) + p_b;
        rho_g_rk=p_rk/(R*T);
        n_rk= N-s*p_rk^.5 *(1-x);
        if s*p_rk^.5 *(1-x)>=N 
            n_rk=0;
        end 
        rho_rk=((n_rk/rho_g_rk)+(1-n_rk)*((x/rho_c)+((1-x)/rho_m)))^-1;    
        Beta_rk=((1/Beta_r)+((N*R*T*p_rk^(-2)-.5*R*s*T*(1-x)*p_rk^(-1.5)+.5*(x/rho_c + (1-x)/ rho_m)*s*(1-x)*p_rk^(-.5))/(N*R*T*p_rk^(-1)-R*s*T*(1-x)*p_rk^(-.5)+(x/rho_c + (1-x)/ rho_m)-(x/rho_c + (1-x)/ rho_m)*N+(x/rho_c + (1-x)/ rho_m)*s*(1-x)*p_rk^.5)^2)/rho_rk)^-1;
        lambda_rk=(rho_rk*S*A^2)/(H*mu);
    K2=(Qi-lambda_rk*delta_p_rk)*(Beta_rk/(rho_rk*V))*dt;
    delta_p(i+1)=delta_p(i)+.5*(K1+K2);
    p(i)=delta_p(i)+p_b;
    Qo(i)=lambda(i)*delta_p(i);
end
time=(0:dt:(N_iters)*dt)';
dp_analytical=(delta_p(1)-Qi./lambda).*exp(-1.*lambda.*Beta.*time./(rho.*V))+Qi.*time./rho;
Qo_analytical= dp_analytical.*lambda;
Vo=(Qo./rho).*10;
Vo_analytical=(Qo_analytical./rho).*10;
semilogx(time, Vo, 'LineWidth', 2)
hold on
semilogx(time, Vo_analytical, '. ')


			Complete