学习中心

Modeling with PDEs: Global ODEs and DAEs Interface

---

In Part 6 of this course on modeling with partial differential equations (PDEs) in COMSOL Multiphysics^®, you will learn how to solve time-dependent ordinary differential equations (ODEs) and differential-algebraic equations (DAEs) that have no space dependency and hence no association with a mesh or geometry. For this purpose, you can use theGlobal ODEs and DAEsinterface. Such equations are also known asinitial value problems. We will also cover linear and nonlinear equations, including systems of equations.

Table of Contents

1. Global ODEs and DAEs
2. Solving an ODE for a Forced Mass-Spring-Damper System
3. Solving a System of ODEs
4. Solving a System of Nonlinear ODEs
5. Plotting Phase Portraits
6. Differential-Algebraic Systems
7. Coupling PDE Models with ODEs
8. The Lorentz Equation System
9. Further Learning

Global ODEs and DAEs

Just as in the case of PDEs, if you have a common application that requires modeling an ODE as an initial value problem, you might find a user-friendly, built-in interface in one of the add-on products. To discover which built-in interfaces are available, start a new model and, in theSelect Space Dimensionwindow, select the0Doption. The term0Dis used for initial value problems that do not have any spatial dependency.

The0Doption.

On the followingSelect Physicspage, you will find a list of available 0D interfaces, as shown in the figure below. To see a complete list of these interfaces, you will need access to various add-on products.

The built-in 0D physics interfaces for modeling ODEs that are initial value problems.

In the figure above, the general-purposeGlobal ODEs and DAEsinterface is highlighted. This is the interface you should select if you want to solve your own ODE or system of ODEs. This interface can also be used if you have a differential-algebraic system of equations, which is a mix of differential and algebraic equations.

Solving an ODE for a Forced Mass-Spring-Damper System

As an easy example to start with, we will look at how we can solve the following second-order initial value problem for a classical forced mass-spring-damper system:

This way of writing a differential equation usesLeibniz's notation.

We can more compactly write this as:

This is known asLagrange's notation.

Yet another option is Newton's notation:

Here, we will use Lagrange's notation.

The forced mass-spring-damper system is an equation in a time-dependent variable,, representing the displacement from an equilibrium position of a point with mass,. The coefficientrepresents damping, andis a spring constant. The right-hand-side functionrepresents an external forcing function. Associated with a second-order ODE are two initial conditions:

For more information on the mass-spring-damper model, seehttps://en.wikipedia.org/wiki/Mass-spring-damper_model.

As a concrete example, we will simulate a mass-spring-damper system with the following parameters:

Assuming you have selected theGlobal ODEs and DAEsinterface, as a final step in the Model Wizard, select aTime Dependentstudy, as shown below.

ATime Dependentstudy is needed for solving an initial value problem.

This will make theGlobal Equationsinterface available, as shown in the figure below.

The mathematics interface for ODEs is calledGlobal Equations.

Before defining the equation, create the following global parameters:

The parameters used for the mass-spring-damper system.

In the user interface, you can use the syntaxutanduttfor the first and second derivatives, respectively.

TheGlobal Equationsinterface requires that we enter the equation on the form:

To adapt the equation to the format required by the user interface, we rewrite the equation slightly as:

so that in all we have:

where the choice of sign here is a convention used in COMSOL^®. The convention is that an external (generalized) force acting on a system should enter the equations with a positive sign.

Note that the interface will also accept nonlinear ODEs on implicit form. As the form of the functionimplies, an ODE can be defined using a mathematical expression that involves the solution variable, its first- and second-order time derivatives, and an explicit dependency on the time variablet.

In theGlobal Equationsinterface, you can now enter the following:

• Name:u
• Equation f(u,ut,utt,t):A*cos(w*t)-m*utt-c*ut-k*u
• Initial value:u(0)=2
• Initial value:u'(0)=0
• Dependent variable quantity, Displacement:m
• Source term quantity, Force:N

In the figure below, these settings have been entered in the user interface:

TheGlobal Equationssettings for the mass-spring-damper ODE.

The next step is to define the time interval that we would like to solve for.

In the Model Builder, select theStudy>Step 1: Time Dependentnode. Change the last time in the expression forOutput timesto 50. Therange(0,0.1,50)syntax means that the solver will output results in intervals of 0.1 s from 0 s to 50 s. These are not the times that the solver is using in its time-stepping algorithm. The algorithm is adaptive and will take more, or sometimes less, time steps than what is specified asOutput timesdepending on the behavior of the equation and the set tolerance.

TheOutput timessettings for solving between 0 and 50 seconds.

ClickComputeto solve. The default plot shows the solution versus time.

The Model Builder with a Global plot selected, the corresponding settings, and the Graphics window showing the mass displacement. The solution to the mass-spring-damper ODE between 0 and 50 seconds.

To get a sense of the level of accuracy of the solution, you can lower the tolerance. In theTime Dependent Settingswindow, change theToleranceoption fromPhysics controlledtoUser controlled, and set theRelative tolerancevalue to 0.0001 (or 1e-4), as shown in the figure below.

TheRelative tolerancesettings for theTime Dependentstudy.

The solution changes slightly, as shown in the figure below.

The Model Builder with a Global plot selected, the corresponding Settings window, and the Graphics window showing the updated mass displacement. The solution to the mass-spring-damper ODE with a lower solver tolerance.

You can even create a parametric sweep over varying tolerance values, as shown in the figure below.

The Model Builder with a Global plot selected, the corresponding settings, and the Graphics window showing different mass displacement results. The solution to the mass-spring-damper ODE with a sweep over successively lower solver tolerances.

Note that the corresponding files are available for download.

We may also be interested in the frequency content of the solution. To display the frequency spectrum in the original mass-spring-damper model (without the tolerance sweep), in theGlobalplotSettingswindow, change theParameteroption in thex-Axis Datasection toDiscrete Fourier transform. Then change theShowoption toFrequency spectrum. This plot is shown in the figure below:

The Model Builder with a Global plot selected, the corresponding Settings window, and the Graphics window showing frequency versus Fourier coefficient. The frequency spectrum of the mass-spring-damper model.

You can plot the axes on a logarithmic (log) scale by clicking thex-Axis Log Scaleandy-Axis Log Scalebuttons in theGraphicstoolbar. This gives the following plot:

The Model Builder with a Global plot selected, the corresponding settings, and the Graphics window, with the graph plotted with a log scale. The frequency spectrum of the mass-spring-damper model, plotted using log scale.

We see that there are two frequency peaks at 0.16 Hz and 0.32 Hz. Note that if you solve for a longer time period, there will be a sharp spike corresponding to the forcing frequency since it determines the steady-state periodic solution.

Solving a System of ODEs

TheGlobal ODEs and DAEsinterface can also be used to solve systems of ODEs. Within chemical engineering it is common to solve systems of ODEs with hundreds, sometimes thousands, of equations. The Chemical Reaction Engineering Module has specialized modeling features for entering such systems. Within electrical engineering you might likewise have large ODE systems when modeling electrical circuits. This is applicable to theElectrical Circuitsinterface, which is available in several of the add-on products, including the AC/DC Module.

Here, we will look at a simple two-variable system — a system-version of the mass-spring-damper equation. You can rewrite the second-order mass-spring-damper equation as two first-order equations by introducing the velocity,, of the point mass as a new degree of freedom:

The resulting system can be written as two first-order equations:

There are several different ways that we can enter this system in the user interface. We can simply add the equation for the velocity in the second row of theGlobal Equationsinterface, as shown in the figure below.

The Model Builder with the Global Equations interface selected and the corresponding Settings window. The system of differential equations that is equivalent to the scalar mass-spring-damper differential equation.

This system will solve correctly; however, if we do it this way, we run into a problem with handling the units. The fact that we now have inconsistent units is indicated by the yellow-orange coloring of the equations. If we want consistent unit handling we instead need to add a secondGlobal Equationsinterface and change theDependent variable quantityandSource term quantityunits tom/s. This is shown in the figure below.

The second differential equation in the velocity variable used for rewriting the second-order mass-spring-damper differential equation as a system of first-order differential equations.

Solving this system of first-order equations will give us the same solution as when we solved it as a scalar second-order equation.

We can also plot the solution in the phase plane where we plot the displacement,u, on one axis and the velocity,v, on the other axis. This is a 2D plot, and to accomplish this we need to first create a 2D dataset. Right-click theResults>Datasetsnode and, underMore 2D Datasets, selectGrid 2D, as shown in the figure below.

TheGrid 2Ddataset used to create a phase plane plot.

We do not need to change the settings for theGrid 2Ddataset since we will plot quantities that have a global variable scope.

Next, right-click theResultsnode and select a2D Plot Group. Right-click this2D Plot Groupand selectPoint Trajectories, which is available underMore Plots. In thePoint Trajectories Settingswindow, enteruandvfor thex-expressionandy-expression, respectively, as shown in the figure below.

ThePoint Trajectoriessettings used for the phase plane plot.

Next, clear thePlot dataset edgescheck box in the2D Plot Groupsettings, as shown below.

ThePlotdataset edges are cleared for thePoint Trajectoriesplot.

We now get the following phase-plane plot of the solution:

The Model Builder with the Point Trajectories plot selected, the corresponding settings showing the plot data, and the Graphics window showing the plot. The phase-plane plot for the mass-spring-damper differential equation.

We can get a better feeling for the behavior of the solution if we solve the ODE up to 200 s, instead of the current 50 s.

Change theOutput timessetting torange(0,0.01,200)and compute. This result in the following solution:

The Model Builder with the Global plot selected, the Settings window with the y-Axis Data section expanded, and the Graphics window showing the mass displacement.

The solution stabilizes to a sinusoidal form determined by the frequency of the forcing function. For c=0.125[kg/s], the corresponding phase-plane plot looks like this:

The Model Builder with Point Trajectories selected, the corresponding settings, and the Graphics window with the phase-plane plot.

The stable solution orbit is reflected by the darker region in the phase-plane plot.

Note that you can also generate a phase-plane plot for the previous scalar mass-spring-damper equation by enteringutinstead ofvin thePoint Trajectoriesplot. However, the variableutwill then be represented by a lower-order approximation thanu, so the phase-plane plot will look smoother if we instead use the system where the dependent variablesuandvwill have the same discretization order.

Solving a System of Nonlinear ODEs

TheGlobal ODEs and DAEsinterface is not limited to solving linear equations, such as the mass-spring-damper system. We can enter virtually any type of nonlinearity in the equations. A classical nonlinear extension of the mass-spring-damper system type of equation is to add a nonlinear gain or damping term, as follows:

Here, we change the sign of the damping coefficient and multiply it by a factor,. This coefficient causes gain ifis small and causes damping ifis large. In this way, the equation has an interesting interplay between damping and gain. This equation has been extensively studied and is called thevan der Pol oscillator. It has commonly been used to model oscillations in vacuum tubes. However, in our example we will model the van der Pol equation as a continuation of the mass-spring-damper system, even though the physical interpretation is not straightforward in this case.

To enter this equation, we can either modify the second-order mass-spring-damper equation or the first-order system of equations. Here, we will modify the first-order system, as follows:

The first of these equations is shown in the figure below, using the same sign convention as in the earlier example.

A close-up of the Model Builder with the Global Equations interface selected and a close-up of the Global Equations table in the Settings window.

The differential equation for the van der Pol oscillator, entered in theGlobal Equationsinterface.

The coefficienteis used to generate a consistent unit and is defined as shown below:

The parameters for the van der Pol oscillator differential equation.

Solving this, usingOutput timesset torange(0,0.01,200), gives the following solution:

The Model Builder with the Global plot selected, the settings, and the Graphics window showing the solution.

The solution to the van der Pol oscillator differential equation.

It also gives the following phase-plane plot:

The Model Builder with the Point Trajectories plot selected, the corresponding Settings window, and the Graphics window showing the van der Pol oscillator phase plot.

The phase plot for the van der Pol oscillator differential equation.

The darker portion of this plot indicates that the solution has a so-calledlimit cycle. For the solution in the time domain, this is due to the fact that it approaches a limiting periodic solution. However, this limiting waveform is not sinusoidal but has a more complicated behavior. We can see its shape more clearly by zooming in on the last few cycles in the time-domain plot, as shown below:

The Model Builder with the Global plot selected, the plot settings, and the Graphics window with the plot zoomed in.

A zoomed-in view of the solution to the van der Pol oscillator differential equation.

By changing the value of the damping (or gain) coefficientc, you can get other waveforms. Note that the nonlinearity of the equation also shifts the dominant frequency. The figure below shows the waveform for c=10[kg/s].

The Model Builder with the Global plot selected, the Global Settings window, and the Graphics window showing the new mass displacement waveform.

The solution to the van der Pol oscillator differential equation for a different damping coefficient.

This waveform has a much more complicated frequency spectrum than the linear mass-spring-damper system, as shown in the figure below using logarithmic scaling:

The Model Builder with the Global plot selected, the corresponding Settings window, and the Graphics window showing the frequency spectrum become more complex.

The frequency spectrum of the solution to the van der Pol oscillator differential equation.

Plotting Phase Portraits

Many text books have illustrations of thephase portraitof the unforced van der Pol oscillator:

We get an expression that we can use to plot the phase portrait by first rearranging:

We will get the phase portrait by identifyinguwith thex-coordinate andvwith they-coordinate. To get the arrow plot we want, we set:

as thex-component

and

as they-component

We will plot this vector field using aGrid 2Ddataset with coordinatesxandyinstead ofuandv. This means that the actual vector-field expression to use will be:

x component:y

y component:(c*(1-e*x^2)*y-k*x)/m

The figure below shows a phase portrait for c=1[kg/s] and a trajectory for the initial condition (u, v) = (-3, -3).

The Model Builder with the Arrow Surface plot selected, the plot settings, and the Graphics window showing the phase portrait of the van der Pol oscillator solution.

A phase portrait of the solution to the van der Pol oscillator.

A close-up view of the phase portrait of the van der Pol oscillator solution.

Close-up of the phase portrait of the solution to the van der Pol oscillator.

The corresponding file is available in the set of downloadable files associated with this article.

Differential-Algebraic Systems

In Part 5 of this course, we described how a time-dependent Joule heating model is a form of a DEA system since one of the equations does not have time derivatives. In the same way, theGlobal ODEs and DAEsinterface can handle DAE systems.

For example, the following system has a time derivative in the first equation, but not in the second:

The figure below shows this DAE system in the user interface:

The differential algebraic equation system entered in theGlobal Equationsinterface.

Note that theInitial valuefor the second equation does not come into play since there is no time derivative.

The figure below shows the solution to this system:

The Model Builder with the Global plot selected, the settings, and the Graphics window showing the state variable u in blue and the state variable v in green.

The solution to the differential algebraic equation system.

The parameters used in this example are shown in the figure below:

The parameters to the differential algebraic equation system.

Coupling PDE Models with ODEs

You can add one or moreGlobal Equationsinterface to an existing physics interface in order to create customized couplings with ODEs. The figure below shows aGlobal Equationsinterface added to anElectric Currentsinterface:

TheGlobal Equationsettings for an algebraic equation linked to the equation of anElectric Currentsinterface.

Coupling field variables from a physics interface with one or more ODEs is illustrated in a number of our blog posts, including:

• "Using Global Equations to Introduce Fully Coupled Goal Seeking"
• "Introducing Goal Seeking into the Segregated Solver"
• "Accelerating Model Convergence with Symbolic Differentiation"

In addition, a related topic that you should explore is state variables. For more information on this topic, see our blog post "How to Use State Variables in COMSOL Multiphysics^®".

There are some important sign considerations when creating customized couplings between equations. In COMSOL Multiphysics, the convention is that an external (generalized) force acting on a system should enter the equations with a positive sign. This means that the mass-spring-damper equation should, when coupled with other equations, preferably be written as:

This sign convention does not matter as long as you are only modeling using ODEs. However, it does matter if you, for example, couple this with a solid using aSolid Mechanicsinterface. Such a model could be used for a mass that is glued to a surface on an elastic object. If you use the sign convention with, you can couple using a simple bidirectional constraint. If you instead use, then you need to switch the sign in the constraint force expression. In the latter case, the symmetry of the coupled system will be broken, and you will end up with an equation system that has a nonsymmetric system matrix that is less efficient to solve.

The Lorentz Equation System

In the early 1960s, the Lorentz equation system was developed to model atmospheric convection. Lorenz discovered that the solution exhibited chaotic behavior for certain parameter values and initial conditions. The solution, when plotted in phase space, resembles the figure eight, as shown below. The Lorenz attractor is a set of chaotic solutions of the Lorenz system as viewed in phase space. A COMSOL Multiphysics implementation of this system, with detailed documentation, can be downloadedhere.

The Lorentz system, which is nonlinear, can be written as follows:

with initial conditions:

Here, the parametersa,b, andcare generally positive scalar numbers. Not all solutions are chaotic, but Lorenz found that the values a = 10, b = 8/3, and c = 28 give a system with a chaotic behavior. The settings below show the Lorentz equation system with the following initial conditions:

• u= −10
• v= −4.45
• w= 35.1

The parameterpertis used to illustrate how the system is sensitive to small changes in the initial conditions.

The Lorentz equation system entered in theGlobal Equationsinterface.

The figure below shows the solution between 0 s and 35 s:

The Model Builder with the Global plot selected, the plot settings, and the Graphics window showing the state variable u in blue, the state variable v in green, and the state variable w in red.

The solution to the Lorentz equation system.

The solution in 3D phase space shows the characteristic shape of the Lorentz attractor.

The Lorentz attractor with the outer section being red, the middle section being yellow, and the inner section being green.

The phase portrait of the Lorentz equation system showing the characteristic shape of the Lorentz attractor.

Among the files available for download with this article, you will also find a version of the van der Pol oscillator with a similar type of phase-space plot in 3D, with time (scaled) as the third axis:

A 3D phase-space plot showing the shape in mostly purple and red, with black near the bottom and yellow near the top.

A 3D phase-space plot of the solution to the van der Pol oscillator equation.

Further Learning

Although the models used throughout this article are available for download, we encourage you to build the models yourself following the guidance provided here. This will help to reinforce how to set up ODEs in the software using theGlobal ODEs and DAEsinterface.