Ordinary Differential Equations with Comsol Multiphysics
Denny Otten1
Department of Mathematics Bielefeld University
33501 Bielefeld Germany
Date: April 17, 2015
1. Introduction and Mathematical Setting
In the following we explain how to solve the following ordinary differential equation u′′(t) =− gR2
(u(t) +R)2, t > t0, u(t0) =u0, u′(t0) =v0,
where the gravitational accelerationg (in sm2), the Earth’s radius R (in m), the initial time t0
(in s), the initial height u0 (in m) and the initial velocity v0 (in ms) are given. We seak for a time-dependent functionu(in m) that descibes the height of the body at timet(in s) measured from the Earth’s surface.
For our simulation we use the following values g= 10m
s2, R= 107m, t0= 0s, u0= 0m, v0= 10m s.
2. Model Wizard
Start Comsol Multiphysics.
To start Comsol Multiphysics 5.0 open theTerminaland enter
• comsol -ckl Model Wizard.
• In theNewwindow, clickModel Wizard.
• In theModel Wizardwindow, click0Din theSelect Space Dimensionmenu.
• In the Select Physics tree, select Mathematics>ODE and DAE interfaces>Global ODEs and DAEs (ge).
• ClickAdd.
1e-mail:dotten@math.uni-bielefeld.de, phone:+49 (0)521 106 4784,
fax: +49 (0)521 106 6498, homepage: http://www.math.uni-bielefeld.de/~dotten/.
1
• ClickStudy.
• In theSelect Studytree, selectPreset Studies>Time Dependent.
• ClickDone.
Unit System.
• In theModel Builderwindow, clickUntitled.mph(root).
• In theSettingswindow for Root, locate theUnit Systemsection.
• From theUnit Systemlist, chooseSI.
Some Advanced Settings.
Hint: In theModel Builderwindow you should click on theShowicon and enable everything that is possible from the menu: Expand Sections (Equation View, Override and Con- tribution,Discretization,Stabilization, Advanced Physics Options, Advanced Study Optionsand Advanced Results Options). Done this, clickExpand Allicon.
3. Ordinary Differential Equation
• In theModel Builderwindow, expand theComponent 1>Global ODEs and DAEs (ge) node, then clickGlobal Equations 1.
• In theSettingswindow for Global Equations, locate the Global Equationssection.
• In the table, enter the following settings:
Name f(u, ut, utt, t)(1) Initial value (u_0) (1) Initial value (u_t0) (1) Description
u utt−F u0 v0
• In theSettingswindow for Global Equations, locate the Unitssection.
• ForDependent variable quantity specify Length (m).
• ForSource term quantity specify Acceleration (m/sˆ2)
4. Parameters and Variables
Parameters.
• In theModel Builderwindow, expand theGlobalnode, right-clickDefinitionsand select Parameters. (Alternatively: On theModeltoolbar, clickParameters.)
• In theSettingswindow for Parameters, locate the Parameters section.
• In the table, enter the following settings:
Name Expression Value Description
g 10 10 gravitational acceleration (inm/s2)
R 10ˆ7 1.0000E7 Earth’s radius (inm)
t0 0 0 initial time (ins)
u0 0 0 initial height (inm)
v0 10 10 initial velocity (inm/s)
Variables.
• In theModel Builderwindow, expand theComponent 1 (comp1)node, then clickGlobal Equations 1. (Alternatively: On theModeltoolbar, clickVariablestheLocal Variables.)
• In theSettingswindow for Variables, locate theVariablessection.
• In the table, enter the following settings:
Name Expression Unit Description
F −g∗Rˆ2/(u+R)ˆ2 1/m2 right hand side of 2nd order ODE Hint: Note that thevariablesmust be chosenlocal, not global.
5. Study Settings and Computation
Study Settings.
• In the Model Builderwindow, expand the Study 1 node, then click Step 1: Time De- pendent.
• In theSettingswindow for Time Dependent, locate the Study Settingssection.
• In theTimestext field, typerange(t0,0.1,2).
Computation.
• On theModeltoolbar, clickCompute.
6. Postprocessing and Graphical Output
The time-height plot.
• In theModel Builderwindow, expand theResultsnode, then click1D Plot Group 1.
• In theSettingswindow for 1D Plot Group, locate theTitlesection.
• From theTitle typelist, chooseManual. In theTitletext field, typeVertical throw and free fall of a body.
• Locate thePlot Settingssection. Markboth check boxes,x-axis label andy-axis label.
• In thex-axis labeltext field, enterTime t (s). In they-axis labeltext field, enterHeight u (m).
• In the Model Builder window, expand the Results>1D Plot Group 1 node, then click Global 1.
• In theSettingswindow for Global, locate theLegendssection.
• Clearthe checkboxShow legends.
• The result is shown in Figure7.1.
The time-velocity plot.
• In theModel Builderwindow, right-clickResultsand select1D Plot Group.
• In theModel Builderwindow, expand theResultsnode, then click1D Plot Group 2.
• From theTitle typelist, chooseManual. In theTitletext field, typeVertical throw and free fall of a body.
• Locate thePlot Settingssection. Markboth check boxes,x-axis label andy-axis label.
• In thex-axis labeltext field, enterTime t (s). In they-axis labeltext field, enterVelocity u_t (m).
• In the Model Builder window, expandResults, right-click 1D Plot Group 2 and select Global.
• In theSettingswindow for Global, locate they-Axis Datasection.
• In the table, enter the following settings:
Expression Unit Description
comp1.ut m/s State Variableu, first time derivative
• In theSettingswindow for Global, locate theLegendssection.
• Enablethe checkboxShow legends.
• The result is shown in Figure7.2.
The time-acceleration plot.
• In theModel Builderwindow, right-clickResultsand select1D Plot Group.
• In theModel Builderwindow, expand theResultsnode, then click1D Plot Group 3.
• From theTitle typelist, chooseManual. In theTitletext field, typeVertical throw and free fall of a body.
• Locate thePlot Settingssection. Markboth check boxes,x-axis label andy-axis label.
• In thex-axis labeltext field, enterTime t (s). In they-axis labeltext field, enterAccel- eration u_tt (m).
• In the Model Builder window, expandResults, right-click 1D Plot Group 3 and select Global.
• In theSettingswindow for Global, locate they-Axis Datasection.
• In the table, enter the following settings:
Expression Unit Description
comp1.utt m/sˆ2 State Variableu, second time derivative
• In theSettingswindow for Global, locate theLegendssection.
• Enablethe checkboxShow legends.
• The result is shown in Figure7.3.
7. Save the Model
Save File.
• SelectFile>Save As....
• Select a desired folder, where the model should be saved, and enterODE.mphas theName for the model.
• ClickOK.
Figure 7.1. Heightu(t)(in m) of the body at timet(in s)
Figure 7.2. Velocityv(t) =u′(t)(in ms) of the body at timet(in s)
Figure 7.3. Accelerationa(t) =u′′(t)(in ms2) of the body at time t(in s)