2.7: Solving ODEs with Mathematica- How to find numerical and analytical solutions to ODEs with Mathematica (2024)

Introduction

Mathematica is an advanced mathematics solution program created by Wolfram Research, Inc. One of the most powerful software packages of its kind, Mathematica is able to easily and conveniently solve complicated mathematical formulae, including differential equations. This article focuses on the modeling of first and higher order Ordinary Differential Equations (ODE) in the following forms:

\[\frac{d y}{d x}=f(x, y) \tag{Basic ODE} \]

\[\frac{d y^{z}}{d^{z} x}=f(x, y) \tag{Higher Order ODE} \]

Like all software, Mathematica uses a specific language in order to perform calculations. The names of all functions must be capitalized- this applies to both mathematical functions (such as Sin and Cos) and built-in functions (such as Plot and DSolve). For this reason it is common for users to write their own functions in minimized letters. This decreases the chance of overriding or redefining a Mathematica function. Square brackets always follow function calls, and the function's parameters are always enclosed in curved brackets. For example, if a user wished to plot sin(x) from x = 0 to x = 1, they would type: Plot[Sin[x],{x,0,1}]. The user must type "Shift"+"Enter" to input the function. Typing only "Enter" will add a line to the formula. . For full PC keyboards, the “Enter” key on the far right is the equivalent of “Shift”+”Enter”. If the user does not want to see the output line, it can be suppressed by typing a semi-colon at the end of the expression. Mathematica also recognizes common mathematical constants such as pi (Pi), Euler's constant (E), and imaginary numbers (I). Note that these constants must also be capitalized when entered.

Mathematica features two functions for solving ODEs: DSolve and NDSolve. DSolve is used when the user wishes to find the general function or functions which solve the differential equation, and NDSolve is used when the user has an initial condition. The prompts for the functions are quite similar. Note: Two equal signs must be used to denote equality in an equation. Using one equal sign assigns a value to a variable.

Example:

f(x) = 5x² + 7 This input creates a variable named f(x) that has a value of 5x² + 7.

f(x) = = 5x² + 7 This input creates a function f(x) which is defined by 5x² + 7.

Mathematica will sometimes store user created functions or variables into its cache. This causes some (potentially) correct commands to output errors. It is recommended to quit the local kernel and reinitialize the commands after any changes to your code.

Also, check that every command in the code is being executed. Clearing all the outputs in the cell may be helpful to figure out which commands had not been executed. This can be done by going to the “Cell” option found on the top and choosing the “Delete All Output” Example:

OUTPUT IS NOT CLEARED

CLEARING THE OUTPUTS

CODE WITH CLEARED OUTPUTS

As seen in the examples, it is easier to troubleshoot and debug your program when it looks less confusing. Clearing the extra outputs helps you focus on just the code that you have written.

First Order ODEs

Notation: In the following examples, eqn represents the ODE, y represents the function being solved for, i represents the initial condition, and x and t are independent variables.

Second Order and Higher ODEs

Input for Higher Order ODEs: Solution = NDSolve[{eqn,y[0] == i1, y' [0] == i2,…},y,{ x,xmin,xmax}]

Plotting the solution: Plot[Evaluate[{y[x], y' [x],...} /. Solution],{x,xmin,xmax}]

Here is a simple example of what you would type into Mathematica:

• Solution = NDSolve[{y[x] +Sin[y[x]]+y[x]==0,y[0]==1,y'[0]==0,y,{x,0,30}]

Mathematica will output:

• Output[1]= {{y->InterpolatingFunction[ { {0.,30.} },<>] } }

To plot this function you type:

• Plot[Evaluate[{y[x],y'[x],y[x]}/.Solution],{x,0,30}]

Mathematica should output a graph.

Algorithms Used by Mathematica

Mathematica uses two main algorithms in order to determine the solution to a differential equation. These algorithms are the Adams method and the Gear method.

The Adams and Gear methods are forms of linear multistep methods. An example of these would be the following:

In the example above, h denotes the step size and the coefficients are determined by the method used. Multistep methods are expansions of more familiar single-step methods used to solve differentials (i.e. Euler, Runge-Kutta, Taylor). Each of these methods requires an initial point in order to calculate the following point. Similarily, multistep methods also require initial points in order to solve the ODE. The number of initial points required depends on which method is used to solve the ODE. Multistep methods typically produce less error than the single-step methods because of multiple initial points.

In order to determine what method to use one must first find the stiffness of the function.

Worked out Example 1

You have been assigned a non-isothermal reactor with the following material and enthalpy balances.

\[\\frac {dc_a}{dt} = c_a(t)^2-c_b(t)- T(t)^2 \nonumber \]Species a Material balance

\[\\frac {dc_b}{dt} = c_b(t)^3 - c_a(t)- T(t)^2 \nonumber \]Species b Material balance

\[\\frac {dT}{dt} = c_a(t) - c_b(t) - T(t)^2 \nonumber \]Enthalpy balance

When is the concentration of species a and species b equal? The initial conditions for the reactor are as follows:

T₀ = 1

c_a0 = 2.0

c_b0 = 1.8

t_f = 0.2

s = NDSolve[{x'[t] == x[t]^2 - y[t] - z[t]^2, y'[t] == y[t]^3 - x[t] - z[t]^2, z'[t] == x[t] - y[t] - z[t]^2, x[0] == 2, y[0] == 1.8, z[0] == 1}, {x, y, z}, {t, 0.2}]

Plot[Evaluate[{x[t],y[t]} /. s], {t, 0, 0.2}]

Worked out Example 2

You have been asked to study a semi-batch reactor for the reaction . The rate law is , where k = 2.7. Other parameters are: initial volume = 5, volumetric flow rate = 0.05, initial concentration of A = 2500. Use Mathematica to create a Conversion versus Time plot for 100 time units.

1. Mole Balance:

\[\\frac{dX}{dt} = - \frac{r_aV}{N_{A0}} \nonumber \]

2. Rate Law:

3. Stoichiometry:

4. Combine:

If the semi-batch reactor runs for 100 time units, the conversion is about 0.8. This can be seen in the plot below or by running the two commands that follow into Mathematica.

s = NDSolve[{y'[t] == 2.7 (1 - y[t]) (5 + .05 t)/(2500 y[t]), y[0] == 0.0001}, y, {t, 0, 100}];

Plot[Evaluate[y[t] /. s], {t, 0, 100}, PlotRange -> All]

This is the Mathematica notebook file for the example: Media:MathematicaEx2.nb

Note: The parameters used in this problem are fabricated and are intended to illustrate the use of Mathematica in solving ODEs.

Formatting Plots in Mathematica

Here are some useful tips to formatting plots in Mathematica:

Label Plot: e.g. Plot[Sin[x], {x,0,2Pi}, PlotLabel->"Volume of Tank vs. Time"]

Label Axes: e.g. Plot[Sin[x], {x,0,2Pi}, AxesLabel->{"time(min)","volume(L)"}]

Color Plot: e.g. Plot[Sin[x], {x,0,2Pi}, PlotStyle->Red]. You can also make the plot thick, dashed, etc. e.g. Plot[Sin[x], {x,0,2Pi}, PlotStyle->{Red,Thick,Dashed}]

Insert Legend: e.g. Needs["PlotLegends`"] <--MUST insert this BEFORE your plot command. Plot[{Sin[x],Cos[x]}, {x,0,2Pi}, PlotLegend->{"sine", "cosine"}]

Sage's Corner

Solving ODEs with Mathematica Video

video.google.com/googleplayer.swf?docId=5007050669435611182

slides for this talk

Additional Tips and Tricks for Troubleshooting in Mathematica

Mathematica is a powerful computing tool, however the syntax can be a bit difficult to understand. Here are some notes for troubleshooting in Mathematica.

1. Check to make sure that your variable names and signs are consistent.

• Ex) Make sure you use xI everywhere instead of xI and x1 or xl.
• Ex) Functions, including the ones you create, are usually followed by brackets such as Sin[x] or y[x]. However, brackets are not necessary when you are solving for a function in a set of differential equations such as NDSolve[eqns, {y}, {x, 0, 50}];
• Ex) Check to see if your parentheses are aligned such that you are actually entering the function you think you're entering. Recall order of operations and the distributive property. x*(a+b) is NOT equal (x*a) + b. This seems simple, but often gets overlooked when dealing with several lines of code.

2. You may find it easier to define all of your variables together at the beginning of your code. This makes it easy to go back and change the values assigned to each variable if trying to determine the impact one variable has on a system of equations. For instance, say you are trying to determine the effects the flow rates into two different tanks (F1, F2) will have on the tank volumes over ten time steps. The differential equations governing the situation are: and , where F1 = 2, F2 = 12, V1(0) = 0, V2(0) = 0.

If you write the ODEs in Mathematica by directly substituting in F1 = 2 and F2 = 12, you will have to change the ODEs each time you change the values of F1 and F2. Below is an example of what the Mathematica equations would look like.

s = NDSolve[{V1’[t] == 8 – V1[t], V2’[t] == 12 – (1/3)*V2[t], V1[0] == 0, V2[0] == 0},{V1,V2},{t,0,10}]

Another option is to define F1 and F2 before your equations, and then call them when solving for your ODEs. This allows you to easily substitute in new values for your variables without changing the differential equations, and reduces the odds of making simple computational errors in your ODE in the process. The Mathematica code would look as shown below.

variables = {F1 -> 2, F2 -> 12};

s = NDSolve[{V1’[t] == 4*F1 – V1[t] , V2’[t] == F2 – (1/3)*V2[t], V1[0] == 0, V2[0] == 0} /. variables,{V1,V2},{t,0,10}]

3. Be Aware of the Kernel

The Mathematica Kernel stores all of the information about variable definitions. When you define a variable , the definition is stored there. The kernel is automatically started when a new Mathematica session is started. You may also start it manually. To start the kernel manually go to Evaluation -> Start Kernal -> Local. Once the kernel is started and you wish to go back and change a variable definition you must "Quit the Kernal" before you see the change occur. In order to "Quit the Kernal" you must go to Evaluation -> Quit Kernal -> Local. Terminating the Mathematica Kernal, erases all of the definitions previously entered. For this reason, after you "Quit the Kernal" and enter in new definitions for your variables you must re-enter all of your code. The images below show how to "Start Kernal" and how to "Quit the Kernel."

Example:

An example of when you would want to "Quit the Kernal" when using Mathematica for Controls is when you are finding a steady state point using a PID controller. We start by defining the variables for vset, taui, taud, and Kc. Mathematica solves this equation using the defined variables and shows that (x->3.9, V->10) is the steady state point. This is shown in the figure 1 below.

Now you see that the definition for Kc has been deleted because the steady state points are in terms of Kc. To find the impact of Kc you can use this solution.

4. Define functions or formulas that are used often

If you will use a function often and it is not already defined in Mathematica by default, define your own. This can be especially helpful if you intend repeat the same function multiple times while only changing variables. A scenario where defining your own function could save you time and error is when comparing P-Fisher's values while keeping the margins constant. Only the variables will change but the function will remain the same. The new function will only last as long as each session in Mathematica.

To define a function, begin by naming the function followed by an open bracket. List all the variables with an underscore after each name and separate each variable with a comma. Close the bracket when all the variables have been listed. To let Mathematica know that this is a newly defined function, after closing the bracket place a semi-colon and equal sign. Now define what the formula is using the variable names. When you have finished, click shift-enter simultaneously. To test the new formula, start with naming the function and open bracket, then list all the variables numerical values separated by comma's in the same order that was listed when defining the function. Close the bracket after all the variables have been listed and click shift-enter simultaneously. The output will be the next line in Mathematica with the answer.

Example:

It is also possible to use this tool to put equations in terms of unknown variables. To do this begin by defining the function the same way as before. When plugging in the numerical values of the variables, leave the unknown variables as their variable name. The Mathematica output will provide the answer in terms of the unknown variables.

Example:

5. "Make it pretty"

Oftentimes when people program in any language, they tend to crowd lines, almost as if they were trying to save paper. Well, since this does not apply when you are typing on a screen, don't worry about it.

• Ex) You could write:

vars = {V -> 1000, Cao -> 0.2, UA -> 20000, F -> 2000, DE1 -> 30000, k1 -> 1.0 10^13, DH1 -> 50000, rcp -> 60, R -> 1.987, To -> 500, Tf -> 500, Tjin -> 500, Fj -> 1000, Vj -> 100, Caf -> 0.2, DH2 -> -192000, DE2 -> 30000, k2 -> 1.0 10^13};

eqns = {Ca'[t] == (Caf - Ca[t]))/V - Ca[t]*k1*Exp[-DE1/(R*T[t])], Cb'[t] == (0 - Cb[t])/V +k1*Ca[t]*Exp[-DE1/(R*T[t])] - k2*Cb[t]*Exp[-DE2/(R*T[t])], T'[t] == (Tf - T[t])/V + (-DH1/(rcp))*k1*Ca[t]*Exp[-DE1/(R*T[t])] + (-DH2/rcp )*k2*Cb[t]*Exp[-DE2/(R*T[t])] - (UA (T[t] - Tj[t]))/(V *rcp), Tj[t] == (Fj (Tjin - Tj[t]))/Vj + (UA (T[t] - Tj[t]))/(Vj *rcp), Ca[0] == 0.2, Cb[0] == 0, T[0] == 500, Tj[0] == 500};

sol = NDSolve[eqns /.vars, {Ca, Cb, T, Tj}, {t, 0, 50}]; Plot[{Ca[t]} /. sol, {t, 0, 50}, AxesLabel -> {t, Ca}, PlotRange -> Full];Plot[{Cb[t]} /. sol, {t, 0, 50}, AxesLabel -> {t, Cb},PlotRange -> Full];Plot[{T[t]} /. sol, {t, 0, 50}, AxesLabel -> {t, T}, PlotRange -> Full];Plot[{Tj[t]} /. sol, {t, 0, 50}, AxesLabel -> {t, Tj}, PlotRange -> Full]

• but it looks much better if you do this:

These thoughts will hopefully help you later on and make it easier to read your code.

6.) "Check the colors!"

Here is a list of font color that Mathmatica will output for a specific reason.

Local Variables in a certain Font Color

Local variables of Module and With in Green

1. Example:

Function arguments and pattern names in Green (Italics)

1. Example:

Variables made special by use in arguments in Turquoise

• Example:

Errors and Warnings in a certain Font Color

Syntax Errors in Purple

• Example:

Emphasized Syntax Errors in Dark Red with Yellow Background

• Example:

Missing arguments in Bright Red

1. Example:

Excess arguments in Bright Red

1. Example:

Possible unwanted assignments in Bright Red

• Example:

Unrecognized option names in Bright Red

• Example:

Local Scope conflicts in Dark Red with Yellow Background

• Example:

Variables that will go out of scope before being used in Bright Red

• Example:

Shadowing in multiple contexts in Bright Red

• Example:

Other in a certain Font Color

Comments in Light Grey

• Example:

Strings in Dark Grey

• Example:

Global symbols that have no value assigned in Bright Blue

• Example:

Accessing Mathematica from your personal computer

If you prefer to use Mathematica from your home computer rather than one on campus, you can remotely log into a CAEN computer through U-M’s Virtual Sites service. Follow the steps outlined below:

2. Visit http://virtualsites.umich.edu/ . The U-M WEBLOGIN page will appear. Log in using your unique name and Kerberos password.
3. The Virtual Sites page should appear on your screen. Click on the large, orange Connect Now button.
4. Under the heading What software do you need?, select Engineering (CAEN) Legacy XP. Then click on the Request Connection button toward the bottom of the page.
5. Follow the instructions under the heading Launch Virtual Sites by downloading and opening the file in the Internet Explorer information bar. A window will pop up on your screen titled Remote Desktop Connection. Click the Connect button.
6. You will then see the screen that appears when logging into a CAEN computer. Your unique name should already be included. Type your Kerberos password and click OK to log in to the CAEN computer.
7. Click OK on the Notice window that pops up and your personal settings will be applied. Also click OK if a System Requirements Wizard window appears.
8. If you scroll down using the scroll bar on the right-hand side of the screen, you will see the Start button. Click Start -> All Programs -> Math and Numerical Methods -> Wolfram Mathematica -> Mathematica 6. Once Mathematica opens on the computer screen, you may use it as though you were actually sitting in front of a CAEN computer on campus.
9. To log out, click Start -> Log Off. Under the Log Off Windows window, click Log Off.

References

• "Adams Method", Wolfram MathWorld, Online: August 5, 2007. Available http://mathworld.wolfram.com/AdamsMethod.html [1]
• "Gear Predictor-Corrector Method", Davidson University, Online August 5, 2007. Available webphysics.davidson.edu/Projects/SuFischer/node47.html [2]
• Fogler, H. Scott (2006), Elements of Chemical Reaction Engineering, New Jersey: Pretince Hall PTR. ISBN 0-13-047394-4 [3]
• "Adams-Bashforth", Pitt Math, Online: September 8, 2007. Available: www.math.pitt.edu/~sussmanm/2071Spring07/lab02/index.html#AdamsBashforthMethods [4]
• "Lipschitz Condition", Wolfram MathWorld, Online: September 7, 2007. Available http://mathworld.wolfram.com/LipschitzCondition.html [5]
• "Adams-Bashforth (AB) Methods", University of California-San Diego, Online: September 7, 2007. Available: renaissance.ucsd.edu/chapters/chap11.pdf [6]
• "Gears Method", University of California-San Diego, Online: September 7, 2007. Available: renaissance.ucsd.edu/chapters/chap11.pdf [7]
• Use the Harvard Referencing style for references in the document.
• For more information on when to reference see the following Wikipedia entry.