The material form Mathematica- Wolfram lctures in Pchem

How to define a function

I want to calculate the pressure p of a gas as a function of the temperature T and molar volume v, by using the Beattie-Bridgeman equation of state

	 p = (RT/(v^2)) (1 -  (c/(vT^3))) (v +B-(bB/v)) -   (A/(v^2))(1-(a/v))  	   (1)
	 
Here R is the gas constant. I will perform the calculation for ethane for which 
the parameters a, b, c A, B are  (CHECK!!!!!!!!!!!!!):
	A=2.
	B=3.1
	a=5.4
	b=5.9
	c=34.7
The units of these parameters are such that if I use v in liter/mole, T in degrees Kelvin and R=0.082 liter atm/mol K, I get the pressure in atmospheres. 
	-we plan to use the function later to make a table with the pressure for the temperatures 100, 200, 300, 400, 500, 600, 700, 800, 900,1000 K and v=20 liter/mol and to plot how the pressure varies  with temperature in the interval from 100 K to 1000K (when v is held constant at 20 liter/mol.

Here is how to define a function p[T,v] that calculates the pressure, as a function of T and v, by using the Beattie- Bridgeman equation Eq. (1).  (use comments * *)

Clear[pBB]
pBB[T_(*K*),v_(*liter/mol*)]:=Module[{p1,p2,p3, R=0.082 (*atm liter/mol K*)},
    p1=1-c/(v T^3);
    p2=v+B-b B/v;
    p3=A*(1-a/v);
     R T  p1 p2/v^2-p3 /v^2  (*atm*) ]

The definition of the function has a number of elements (shown in different colors) and we will examine them gradually.  In the first line I clear the symbol pBB. It is prudent to do this whenever you define a function, to avoid conflict with possible prior definitions. However, note that if you defined a quantity pBB earlier in the program Clear[pBB] will wipe out all remembrance of this definition. If you clear that definition you must make sure that you no longer need it.

The part in red is the program that calculates the pressure according to the Beattie- Bridgeman formula. Since this expression is fairly long, I break up the calculation into smaller parts. This is a good practice: the program is easier to read and to modify. I give names to various parts of the expression 

	 p = (RT/(v^2)) (1 -  (c/(vT^3))) (v +B-(bB/v)) -   (A/(v^2))(1-(a/v))  
as follows
	 p1=1-c/(v T^3); 
	 p2=v+B-b B/v;
	 p3=A*(1-a/v);
Next I use these symbols to calculate the pressure p, which is given by:

	(R T p1 p2-p3)/v^2 
	
Note that these expressions are placed inside a Module with the following syntax:

	pBB[T_,v_]:=  Module[{expressions },  
	p1=1-c/(v T^3); 
	 p2=v+B-b B/v;
	 p3=A*(1-a/v);
	(R T p1 p2-p3)/v^2 ]
	
	The expressions in red constitute the program (sometimes called the body of the function) that calculate the pressure according to Eq. (1).  In the first red line the program calculates p1 and gives no output (because I used ; after the line).  It then goes on and calculates p2, then p3, then (R T p1 p2-p3)/v^2. The last expression is the output (it does not have ; after it, so the computer prints on the screen the result of evaluating it). 
	In the program (the red lines inside the Module) you must use ; after every line of code (at the end of every expression) except the last (which is the output).  If you forget ; you are lucky if the computer gets confused and gives you an error message. Otherwise it is possible that it will not calculate what you intended.

The stuff in blue in the expression below

	pBB[T_,v_]:=  Module[{p1,p2,p3,R=0.082 },  
	p1=1-c/(v T^3); 
	p2=v+B-b B/v;
	p3=A*(1-a/v);
	(R T p1 p2-p3)/v^2 ]
	
tells the computer that the quantities p1, p2, p3 need not be communicated to the program in which the function is used. The assignments are made and remembered only when the red lines are calculated. After that, they are wiped out from the memory (the procedure is more clever, but the outcome is as if the symbols were cleared). The same is true for the assignment R=0.082. This means that the assignments to p1, p2, p3 will never interfere with the program in which the function is used. By placing these quantities in the curly brackets in Module[{ }, ] we have limited their scope to the red lines in the Module. This is computerese for saying that no other part of the program, other then the red lines, will ever use these assignments. Thus we ensure that quantities defined in our function do not "pollute" the program where the function is used.

Note that in several places in the program 

pBB[T_(*K*),v_(*liter/mol*)]:=  Module[{p1,p2,p3,R=0.082 (*liter atm/mol K*)},  
	p1=1-c/(v T^3); 
	p2=v+B-b B/v;
	p3=A*(1-a/v);
	(R T p1 p2-p3)/v^2(*atm*) ]
	
I have statements of the form:

	(* something*)
	
The syntax (* something*) is used for comments. The computer ignores them (does not try to interpret them as programs).  A good programmer puts many comments in the program, so that he, or someone else, can read it and use it easily a year after it was written. In my case I used the comments to remind me the units of various quantities.

The left hand side 

	pBB[T_(*K*),v_(*liter/mol*)]= right-hand-side
	
is very interesting.  It says that the expressions right-hand-side define a function named pBB.  T and v are function arguments. You can define a function with one, or two, or three... arguments, depending on what you need. The arguments are always enclosed in square brackets and are separated by commas.
	The syntax 
	
		pBB[T_,v_] 
		
may puzzle you. The symbol  _   is called by typesetters an underscore. In Mathematica it is called a blank.   In the definition of the function  every argument must have an underscore at the right of its name (i.e. T_ not T, v_, not v). It is a very bad error, that could cause a lot of trouble later, to use p[T,v]:= right-hand-side  instead of p[T_,v_]:=  right-hand-side.  Do not forget the underscore!!!  
Here is what the underscore does. The function is

pBB[T_(*K*),v_(*liter/mol*)]:=  Module[{p1,p2,p3,R=0.082 (*liter atm/mol K*)},  
	p1=1-c/(v T^3); 
	p2=v+B-b B/v;
	p3=A*(1-a/v);
	(R T p1 p2-p3)/v^2(*atm*) ]
	
If I type pBB[T,v] anywhere in a program, the computer will evaluate it and replace the symbol pBB[T,v] with expression obtained by evaluating the right-hand-side of the definition of the function. See below how this works:

pBB[T,v]

?pBB

What happens if I want to calculate the pressure for T=300 and p=10?  No problem. Use:

pBB[300,10]

Here is how this works. The computer "reads" the first argument (300 in our example) and replaces T in the right-hand-side of the definition with 300. Then reads the second argument and replaces p in the right-hand-side of the definition with 10. Then it evaluates the program with these values for T and p. The computer does that because of the underscore in pBB[T_,p_].   This tells the machine the following.  If the left-hand-side in the definition of the function is pBB[T_, p_], and you encounter the expression pBB[expr1,expr2], than  replace T in the right-hand-side of the definition with  expr1, and replace p in the right-hand-side with expr2.  Let us see how this works:

pBB[expr1,expr2]

pBB[Sin[x],Tan[y]]

If I use p[300,20] the function replaces T with 300 and v with 20 in the right hand side and evaluates it.

pBB[300,20]

To emphasize the importance of the underscore let us define two simple functions 

	f[x]:=x^2 
and 
	ff[x_]:=x^2. 
	
One does not use an underscore and the other does.

f[x]:=x^2
ff[x_]:=x^2

Let us see if they work correctly:

{f[x],ff[x]}

They both seem to be doing fine! But look at the next evaluation:

{f[2],ff[2]}

The function that did not use underscore works only when the argument is x, but not when the argument is 2, or 3, or any other number or expression. f[x] is useless.

Finally, note that I did not use the equal sign for the assignment. I used instead 

	p[T_,v_]:= something,
	
 where := is called a delayed assignment. To explain what this means I give a simple example. Define the function 
 
	g[x_]:=Sin[x]^3

After typing this I press Shift-Enter and nothing happens (see below)

g[x_]:=Sin[x]^3

An input was made (the In[number] appears in front of the definition) but there is no output. This happens because := tells the computer to memorize the definition of the function but not evaluate it. How do I evaluate the function? The function is evaluated whenever it is encountered in a program. For example, if I type g[3] and the press Shift-Enter I obtain

g[3]

This is why := is called a delayed assignment: it is evaluated later, when needed, not when the definition is communicated to the computer. There are some subtleties about the use of the delayed assignment that will be discussed later.

To use the function in numerical calculations for ethane, I have to give appropriate values to a, b, c, A, B and R. I use a replacement rule.

Clear[A,B,a,b, c, R]
parameters={A->2.,B->3.1,a->5.4,b->5.9,c->34.7 }

I gave the rule the name parameters (you could use any name you want, as long as it has not been used before).

Note that I have cleared the symbols for the parameters before defining the rule. This is very important. If I made a previous assignment A=3.7, the rule becomes {3.7->2., ...etc} instead of {A->2.,....etc}. This is a catastrophe and you are lucky if the computer catches the error. In most cases it will not and it will just give you an erroneous result.

Let us use the function and the rule to calculate the pressure for ethane:

pBB[T,v]/.parameters

If I want to calculate the pressure at T=200 and v=20, for ethane, I use.

pBB[200,20]/.parameters

At 200 degrees the gas having a molar volume of 20 liter/mol has a pressure of 0.90595 atm.

If I want to use the function pBB to calculate the pressure for other molecules, I have to give it the parameters for those molecules. I do that by changing the values of A,B, a, b, c in the replacement rule

If you want to execute a set of operations many times, you should always make a function that executes those operations. In our case all the operations required by the Beattie-Bridgeman equation were "wrapped" under one symbol, pBB[T,v].This is easy to remember and you do not have to repeat the entire code whenever you need to calculate the pressure.You just write pBB[T,v] and the machine finds the code and uses it.

Creating functions  is a very important part of code writing, so let us practice it a bit. I'll write a function that calculates the pressure according to the Benedict-Webb-Rubin equation of state

	p=RTρ + (BRT-A-(Subscript[C, 0])(T^(-2)))(ρ^2) + (bRT-a)(ρ^3) + a α (ρ^6) +c (T^(-2)) (ρ^3) (1 +γ (ρ^2)) exp[-((γρ)^2)]
	
 where 
	ρ=1/v 
	
and v is the molar volume, and A, B, (Subscript[C, 0]), R, a, b, α and γ are parameters whose value change from molecule to molecule. Here is my proposal for writing this function:

pBWR[T_,v_]:= Module[{ρ,p1,p2,p3,p4,p5,R=0.082 (*liter atm/mol K*)},
Clear[A,B,C0,a,b,α,γ];
 ρ=1/v;
p1=R T ρ;
p2=(B R T-A-C0 T^-2)ρ^2;
p3=(b R T-a)ρ^3;
p4=a α ρ^6;
p5= c T^-2 ρ^3 (1+γ ρ^2) Exp[-γ ρ^2];
 p1+p2+p3+p4+p5 ]

The name of the function is pBWR. The arguments are the temperature T and the molar volume v.

Here is how I use the function:

pBWR[T,v]

The parameters for methane are:

Clear[A,a,B,b,c,C0,α,γ,R];
paraRule={A->2.4,a->3.1,B->0.9,b->5.2,c->20.3,C0->2.2,α->1.1,γ->0.3,R->0.082}

If the temperature is 400K and the molar volume is 4.2 then the pressure of methane is:

pBWR[400,4.2]/.paraRule

Exercise:  Write a function that calculates the pressure according to the van der Waals equation. 
Exercise: Think of any function in Physical Chemistry and write a function that calculates it.

Cell 2: How to plot a function

Define a function:

Clear[f]
f[x_]:=x^2 Sin[x]

If you want to plot it use:

Plot[f[x],{x,0,2 Pi}]

This syntax 

	Plot[f[x],{x,0,Pi}]), 

indicates that the expression to be plotted is f[x]; the expression {x, 0, Pi} tells the computer that x should vary  from 0 to Pi when the plot is made.  
	 The function must be such that f[x] is a number when x is a number.  Look what happens when this requirement is not fulfilled

Clear[a]
Plot[a x^2,{x,0,2}]

You can put labels on the axes and label the graph.

Plot[f[x],{x,0,Pi},
PlotLabel ->"My function",
AxesLabel->{"x","f[x]"}]

Sometimes Mathematica cuts off a graph to show in more detail one part of it.

Plot[Exp[-3 x^2] Sin[100 x]/(0.0001 x),{x,-1,1}]

Here the computer cut off the top of the graph to show the wiggly part in more detail. You can force it to show all:

Plot[Exp[-3 x^2] Sin[100 x]/(0.0001 x),{x,-1,1},
PlotRange->All]

Sometimes you want to focus on a part of the graph. For example, you may want to see in more detail what the function looks in a certain range.  Then you can use the PlotRange option, as shown below.

Plot[Exp[-3 x^2] Sin[100 x]/(0.0001 x),{x,-1,1},
PlotRange->{-1,1}];

This shows only that part of the graph for which f[x] is between -1 and +1. Smaller and larger values of the function are not shown.

Exercise: Plot the expression 
	3. (2+ (e^(a ))) and show only those parts of the plot where the function is between zero and ten.  Let  a vary  between zero and 8. Label the graph and the axes

If you have used Cell 1 and the function pBB is still in memory you can plot how the pressure varies with the molar volume at a constant temperature of 200 K for methane

?pBB

Let's see if  the rule giving the parameters is still in the memory (if I isn't, run Cell 1 again)

?parameters

Plot[pBB[200, v]/.parameters,{v,10,40},PlotRange->All,PlotLabel->"         pressure from Beattie- Bridgeman", AxesLabel->{"v liter/mol","p atm"}]

Cell 3: Lists: how to store data and how to act on it.

All computer languages have a special syntax for storing data. Mathematica uses Lists. Here is the syntax that creates a list:

	aList={ expr1, expr2, expr3, expr4}
	
The list consists of curly brackets (braces) containing the list members expr1, expr2, ....  separated by commas. The members could be anything: numbers, expressions, programs or lists. The aList is the "name" of this particular list. You can give a list any name you want. 
	Here are a few examples of lists:

myLis1={1.2, 3.4, 8.452, Sin[x], Log[Sqrt[x]],{a,b,c}}

The first three elements are numbers, the next two are functions and the last one is the list {a,b,c}

myLis2={{a,b,c},{1,2,3}}

That is a list whose elements are two lists.

You can access the elements of a list by using the following syntax
	myList[[n]] 
which returns the n-th element of the list:

myLis1
myLis1[[3]]

The first output shows the list myLis1 and the second shows the third element of the list.

myLis1
myLis1[[6]]

The sixth element of myLis1 is the list {a,b,c}.  If I want to access the second element of the list {a,b,c} I do:

myLis1
myLis1[[6,2]]

As I told you, a list is used to pool together data. For example, if you want to calculate the pressure at v=5.6 liter/mol and  five  temperatures T=  300, 400, 500, 600, 700 I will put the temperatures in one list

temp={300, 400, 500, 600,700}

To calculate the pressure I use the function pBB defined in cell 1.  The function is

Clear[pBB]
pBB[T_(*K*),v_(*liter/mol*)]:=Module[{p1,p2,p3, R=0.082 (*atm liter/mol K*)},
    p1=1-c/(v T^3);
    p2=v+B-b B/v;
    p3=A*(1-a/v);
     R T  p1 p2/v^2-p3 /v^2  (*atm*) ]

The parameters for ethane are

Clear[A,B,a,b, c, R]
parameters={A->2.,B->3.1,a->5.4,b->5.9,c->34.7 }

Here is how the function is used to calculate pressure:

pBB[400,5.6]/.parameters

If I want to calculate the pressure at all the temperatures in the list temp I do:

temp
pressure=pBB[temp,5.6]/.parameters

The result is a list of pressures corresponding to the temperatures in the list temp.  The first pressure, 5.73319 atm is the pressure at a temperature of 400 K (first element of the list temp), etc. 
	As you can see, when the argument of the function pBB is a list, the function acts on each element of the list and returns a list of results. This is not a general behavior. It happens in this case because pBWR is a simple function. In general we have to force the function to act on each element of the list. This is done with a function called Map. Here is the syntax:
	
	Map[aFunction, aList]
	
An example: 

	Map[f,{a,b,c}]

Clear[f]
f[{a,b,c}]
Map[f,{a,b,c}]

When it encounter the above expression, the computer assumes automatically that f is a function with one argument  and applies it to each element of the list. Since the function has not been defined, the output is {f[a], f[b], f[c]}. I can define f and then act again on the list

f[x_]:=x^2 Sin[ 3 x]
Map[f,{a,b,c,d}]

There is a short hand notation for Map that is very extensively used:

f/@{a,b,c,d}

The two symbols in red mean: map the function at the left of /@ into the list at the right of /@

Cell 4: Simple "nameless functions" or "pure functions".

It happens quite often that I need to evaluate a function only in one place in a program. If I have ten such functions, my program will contain ten names that I need to remember. It is easy to get confused and tired and the program is harder to read and understand, a month after I wrote it. Because of this, Mathematica has "nameless functions" that do the job, do not require you to remember their name, can be read where they are used (so program readability is improved ), and are evaluated faster than ordinary functions.

Here is an example of a nameless function that calculates the ideal gas pressure

Clear[w]
w=R #1/#2&

Here #1 and #2 are place holders for the arguments of the function and & at the end tells the computer that everything at the left of & is the definition of a function.

See below how this is used

w[T,v]

It works like any good old function.  However, in doing what I just did, I defeat the purpose of the nameless function, by giving it the name w.  I did not have to do that, I could have done:

R #1/#2&[T,v]

Here are a few additional examples of pure or nameless functions:

Here is a function that "cubes" its argument:

#^3&[x]

The function cubes any argument (i.e. Sin[x] in the example below)

#^3&[Sin[x]]

Here is another one:

Cos[#^2+Log[#]]&[x]

The next one is more complicated: it has two arguments and creates a list:

{#1^2,#2^3,#1-#2}&[a,b]

Pure functions are used very often. One of the most common uses are in statements in which a function is applied to members of a List. This is done with the functions Map and MapThread which we study in the next two cells.

Cell 5: The use of Map to apply functions to each member of a list

In Mathematica data is stored in Lists.  To act on the data we need a syntax that allows a function to act on each element of a list. This is done with Map. In most cases we want to use a nameless function so that we avoid memorizing names of functions and read the program easier.

For example, if I want to calculate the effect of the function Sin[x^2] on each member of a data list I do this

dataList={a,b,c,d,25.};
Sin[#^2]&/@dataList

Here is a more complicated example:

Sqrt[Sin[#]]+Log[#^2]&/@dataList

For example, I may want to calculate the pressure of a gas at the temperatures 300, 400, ...,1000K, and the volume v=22 liter/mole. Let us first do this calculation by using the ideal gas law.  To do this I make a list with the temperatures of interest

tList={300,400,500,600,700,800,900,1000}

If I want to use the ideal gas law

	 p=R T/v, 
	 
with R=0.082 liter atm/mol K and v=22 liter/mole I use

0.082  #/22&/@tList

The Map function applies the nameless function 

	0.082 #/22&  
	
 (0.082  is R in liter atm/mol K and v=22 liter/mol) to every temperature in the list.  The first number in the result above (i.e. 1.11818) is the pressure at the first temperature in tList (i.e. 300 K), etc. Map allows me to do all the calculations in one shot. I don't have to remember the name of the function, and I can see the definition of the function at the point I am using it (improved readability).

Let's take another example. Calculate the pressure by using the van der Waals equation for CO at the temperatures given in the list tList, for v=10 liter/atm. The van der Waals equation for pressure is 

	p= (RT/(v-b)) -(a/(v^2))
	
For CO   a = 1.505 (liter^2)bar/(mol^2) and b=0.03985 liter/mol.  Use v in liter/mol and T in degrees K to get the pressure in bar.

I already defined the list of temperatures in tList

tList

The calculation is:

Clear[a,b,R]
Map[R #/(10-b)-a/10^2&,tList]

To get numerical values I need to replace R, a and b with their values. I do this by using a replacement rule {R->0.083, a->1.505, b->0.093}

Clear[a,b,R]
Map[R #/(10-b)-a/10^2&,tList]/.{R->0.083,a->1.505,b->0.093}

Or you can use:

(R #/(10-b)-a/10^2&/@tList )/.{R->0.083,a->1.505,b->0.093}

Once you get used to /@,  expressions like the one written above are easy to read. Start from /@. Everything at the left of /@ is a function. At the right of /@ is a list into which the function is Mapped.

Exercises: 
(1) Would the syntax 
	Map[(R #/(v-b)-a/v^2)/.{R->0.083,a->1.505,b->0.093})&,tList]
work?
(2) Write a nameless function that takes the logarithm of each element in tList and returns a List with the results. 
(3) Write a nameless function for
	f(x)=x^2 Sin[a x]
     and use it to calculate a list of results for a=2 and x=1,2, ...10

Now let us calculate the pressure for v=20 liters per mole and the temperatures 100, 200, ...,1000 K. I start by making a list of temperatures.  Mathematica provides a function called Table which allows us to make lists. The syntax is 
	Table[ some function of i, {i, imin, imax, incr}]
See below how this works:

Clear [f]
Table[f[i],{i,2,12,2}]

%/.f->Sin

%//N

Here are other examples of using Table:

Table[x,{x,1,10,2}]

Table[x,{x,1,10,Log[1.2]//N}]

You can make a list of arguments and the Sin of the arguments (who needs trigonometric tables anymore?)

Table[{z ,Sin[z]//N},{z,-4,6,0.2}]

Now let us get back to our pressure calculation.  Use Table to make a list of temperatures called te

te=Table[T,{T,100,1000,100}]

To calculate the pressure I use the function p[T,v] and the parameters defined by the replacement rules called parameters (these were defined in Cell 1).  Let's see what they are:

?pBB

?parameters

Now I can calculate the pressure for each temperature in the list te and for v=20 liter/mol, by using the list as an argument

pressures=pBB[#,20]/.parameters&/@te

I turned the function of two variables p[T,v] into a function of one variable

	p[#,20]/.parameters&
	
Now I can use Map ( Map only works with functions of one variable). The result is the pressure at all the temperatures in the list te.

What if I want to vary the volume and the temperature? No problem. Let's say that the temperatures are in the list te and you want to see the pressure at v=10, 15 and 20 liter/mol at each of the temperatures in the list te. Here is how you do that: you take the expression used above for one volume and put it inside a Table to do the calculations for several volumes. Note that the function is in blue and the Table is magenta colored. Map is in red.

res=Table[p[#,v]/.parameters&/@te,{v,10,20,5}]

The output is a list of lists. The first list corresponds to the pressures for v = 10 and all the temperatures in the list te.  For example, p[100, 10] is 0.915019  and p[100, 15] is 0.609515 and p[300, 15] is 1.83993.

Cell 6: The use of MapThread to act on the elements of lists with a function of many variables.

Often we need to act on the elements of several  lists with a function of several variables. We use for this the function MapThread. Here is how this function works

Clear[f];
MapThread[f,{{a,b,c},{A,B,C}}]

The function MapThread

	MapThread[f,{list1, list2}]
 
 has two arguments. The first argument is a function (f in our example). The second argument is a list {list1, list2}that contains in it two lists, list1 and list2. Since the second argument consists of two lists Mathematica assumes that f is a function of two variables. The two lists must have the same number of elements. See what happens if they don't:

MapThread[f,{{A,B},{1,2,3,4}}]

Mathematica complains because the list {A,B} does not have the same length as {1,2,3,4}.

MapThread can be used with more than two lists

MapThread[f,{{A,B},{1,2},{a,b}}]

Since the second argument contains three lists Mathematica assumes automatically that f is a function of three variables.

As you learn to program, you'll find many uses for MapThread.  I'll give you a few examples here. You'll see many more as we go along.

Example 1:  I have two lists, one with experimental values of a quantity and one with the values calculated with some theory. I want to calculate the percent error made by theory.

expValues=Table[n,{n,1,10}]

To generate the "theoretical results" I use the function Random[Real,{0.01,0.03}] which generates random real numbers larger or equal to 0.01 and smaller or equal to 0.03.

calcValues=Table[n+Random[Real,{0.01,0.03}],{n,1,10}]

The error is 

	errorInPercent = (experimental value-theoretical value) 100/experimental value. 
	
Here is how we use MapThread to calculate the error  in percent:

errorPerCent=MapThread[( #1-#2)*100/#1&,{expValues,calcValues}]

I used in MapThread a nameless function 

	(#1-#2)*100/#1&
	
This takes an element from the first list in position #1, and an element from the second list in position #2, and evaluates the expression (list1Element-list2Element)*100/list1Element.  MapThread does this operation for each pair of elements in the list {expValues, calcValues}.

Example 2:  Now I want to make a table with these quantities. I can use a function called TableForm.
Let us see first the two lists:

expValues

calcValues

errorPerCent

TableForm[{expValues, calcValues, errorPerCent}]

This looks awful and if I want to print it will not fit on the paper.  The structure of the list given to the TableForm is not convenient. Here it is what it looks like

{expValues, calcValues, errorPerCent}

TableForm places the first list in the first row of the table, the second list in the second row etc.  See how this works below:

TableForm[{{a,b,c,d,e,f},{A,B,C,D,E,F},{1,2,3,4,5,6}}]

To get a "vertical table" I should have the data organized as in the example shown below:

TableForm[{{a,A,1},{b,B,2},{c,C,3},{d,D,4},{e,E,5},{f,F,6}}]

How do I get from the list

{expValues, calcValues, errorPerCent}

to a list that has the structure

	{{1, 1.02932, -2.931},{2, 2.02472,-1.230609},....{10,10.0116,-0.11599}}  ?
	
Well, we use our trustworthy MapThread:

toTable=MapThread[{#1,#2,#3}&,{expValues, calcValues, errorPerCent}]

I used above the nameless function {#1, #2, #3}&.  This places an element from the first list in the position #1, one from the second list in position #2, and one from the third list in position #3. The MapThread applies this function to all elements in the lists, and produces the list we want. Now TableForm gives a nice vertical table:

TableForm[toTable,TableHeadings->{None,{"ExpValues","calcValues","error%"}}]

Now the table fits in the page. I also used the option TableHeadings to put headings on the table.  You'll learn more about TableForm when we discuss making tables in detail.

Example 3:  Here is another example of use of MapThread. I want to take each element from the list expValues, and subtract the corresponding element from calcValues and then square the difference.  Then I want to add all these squares.  I need to do this when I apply the least squares method to fit data to an equation.

Plus@@MapThread[(#1-#2)^2&,{expValues,calcValues}]

{a,b,c}//FullForm

Here is how this works:

aList=MapThread[(#1-#2)^2&,{{a,b,c},{A,B,C}}]

The function of two variables 

	(#1-#2)^2& 
	
was applied by MapThread first to the pair of variables a and A, next to the pair  b and B, next to the pair c and C.  The function will take a from the first list and subtract A from the second list and square the result, etc. I get a list of the squared differences. I need to add the elements of this list and this is done by Plus@@aList

Plus@@aList

This gives me what I wanted. I can combine the two statements into one:

Plus@@MapThread[(#1-#2)^2&,{{a,b,c},{A,B,C}}]

Example 4: Here is a final example. I want to plot the dependence of the rate constant k on temperature T. Theory (Arrhenius formula) tells me to expect that a plot of  ln(k) versus 1/T is a straight line. I have experimental data for k for various values of T. How do I make this plot? To do this you will have to learn how to plot a list. Then you will have to learn how to process the lists containing data to plot ln(k) versus 1/T.

Data is organized in lists and lists are plotted with the function ListPlot. Here is how ListPlot works. I want to plot a quantity y that depends on a variable x. I fix a value of x and I measure y. Then I take another x and I measure y, etc. Here are the results

x={0.1,0.8,1.5,2.7,3.4,6.3,8.1,9.2}

The results for y are

y={10,12.4,14.6,18.3,19.2,20.9,22.4,23.1}

The normal way of organizing the data in Mathematica is this: the first value in the list y corresponds to the first value of x in the list x, the second value in the list y corresponds to the second value of x in the list x etc.

To plot y versus x with ListPlot I must organize the data as follows

	{{x1, y1},{x2, y2},....{x8, y8}} 

where x1 is the value of x in the first measurement that yielded y1, etc. I can take the lists x and y and reorganize them by using MapThread

toPlot=MapThread[{#1,#2}&,{x,y}]

ListPlot[toPlot,PlotStyle->PointSize[0.02],AxesLabel->{"x","y"}]

The statement

	PlotStyle->PointSize[0.02]
	
makes the dots in the graph a little larger than the default value (which is too small). 

	AxesLabel->{"x","y"}
	
puts labels on the axes. You'll learn more about making pretty plots in future lessons.  
Now let us get back to our rate constant plots.  The data gives the rate constant k

k={1.24, 1.32, 1.81,2.08,2.29,2.76}*10^8 (*mol/(dm^3)sec*);

Make a list with the temperatures at which k was determined

T={292,296,321,333,343,363};

To make the plot I want I need to have a list of the form
	{{Log[1.24*10^8], 1/292},{Log[1.32*10^8], 296}, etc  }

I get this with MapThread

toPlot=MapThread[{100/#2,Log[#1]}&,{k,T}]//N

You can see how this works. The function takes an element from the first list and calculates the logarithm and the element from the second list and calculates its inverse. The results are member of a two member list. MapThread applies this function to all pairs of values from the two lists.

ListPlot[toPlot,PlotStyle->PointSize[0.03],AxesLabel->{"100/T","ln[k]"}]

Review:

Nameless functions: 
	#1^2 +#2^2&
	#1^2  Sin[#2^2]&

Map

Map[f,{a,b,c}]

f/@{a,b,c}

(1/#)&/@{a,b,c}

f[x_,y_]:=x^2 Sin[y];
f[3,#]&/@{a,b,c}

MapThread

MapThread[g,{{a,b,c},{A,B,C}}]

MapThread[f,{{a,b,c},{A,B,C}}]

Part II: Symbolic algebra and calculus with Mathematica

Cell 7:   How to calculate the derivatives of an expression or a function

The derivative of an expression  expr with respect to x is:

	D[expr, x]

The function D has two arguments.The first is an expression (or a function) and the second indicates the variable with respect to which we take the derivative.

Here are some examples:

D[x^2, x]

D[Cos[x],x]

Mathematica can calculate easily the derivative of a function. For example, let us calculate the derivative of pressure with the molar volume, for the function pBB[T,v] defined in Cell 1 (this gives the pressure according to the Beattie- Bridgeman equation). Here is the function:

Clear[pBB]
pBB[T_(*K*),v_(*liter/mol*)]:=Module[{p1,p2,p3, R=0.082 (*atm liter/mol K*)},
    p1=1-c/(v T^3);
    p2=v+B-b B/v;
    p3=A*(1-a/v);
     R T  p1 p2/v^2-p3 /v^2  (*atm*) ]

Call it, to see if it works:

Clear[T,v]
pBB[T,v]

The derivative of this function with v is:

D[pBB[T,v],v]

If you want to get a simplified expression for the output, you can use Simplify or FullSimplify

D[pBB[T,v],v] //Simplify

If you are not happy with the result you may try FullSimplify

D[pBB[T,v],v] //FullSimplify

In this case the two functions (Simplify and FullSimplify) give the same results. Sometimes they do not.

You can calculate the higher derivatives as follows:

D[x^3,{x,2}]

D[Sin[x],{x,3}]

You can calculate the second derivative of p[T,v] with v:

D[pBB[T,v], {v,2}]//FullSimplify

The derivative of p with T is:

D[pBB[T,v],T]

The mixed, partial derivative 
	(((∂^2)p)/(∂T∂v))
is:

D[pBB[T,v],v,T]

Exercise:  Write a function that gives the pressure according to the van der Waals equation. 

	p= (RT/((v-b) )) -(a/(v^2))
	
	Calculate 
	
	κ(T,v)=v(((∂p)/(∂v)))  and  α[T,v]= T(((∂p)/(∂T)))
	
Plot α[T,20] versus T and κ[300,v] versus v. Use a=5.536, b=0.03049 and R=0.083
Calculate α[T,20] for T={100, 200, ... 1000}K. Use a=5.536, b=0.03049 and R=0.083

Cell 8:   How to calculate the integral of a function or an expression.

It is easy to integrate an expression or a function as you can see below

Integrate[x^2,x]

We can check whether the result is correct by taking the derivative of the result

D[%,x]

In the above expression % means the result of the previous evaluation.

Here is another example:

Integrate[Cos[x]^2,x]

D[%,x]

This does not look like Cos[x]^2. However, if we tell the machine to simplify the result (by using //FullSimplify) it puts it in the equivalent form Cos[x]^2

(1/2)+(1/2) Cos[2 x]//FullSimplify

Here is how you check whether the two expressions are identical

((1/2)+(1/2) Cos[2 x]//FullSimplify)===Cos[x]^2

The syntax for an indefinite integral (the limits are not specified) is therefore:

	Integrate[expression, integration variable]

You can calculate definite integrals. The syntax is

	Integrate[expression, {integrationVariable, lowerIntegrationLimit, UpperIntegrationLimit}]

Clear[x,a,b]
Integrate[x^2,{x,a,b}]

Integrate[Cos[x]^2,{x,0,Pi}]

Above, Pi is the symbol for π.

Pi//N

N, in Pi//N, forces the computer to evaluate Pi numerically.

The integration limits do not have to be numbers. They can be symbols or expressions as you see in the two examples given below:

Clear[a,y]
Integrate[x Cos[x],{x,0,y}]

Clear[a,b,x]
Integrate[x Cos[x],{x,Sin[a],Log[b]}]

Mathematica 5 made integration more complicated than it should be. In some cases the function Integrate gets overly cautious and imagines all kind of possible complications. In this case it gives the result and tells you that it is correct if certain conditions are satisfied. Often it is best to ask integrate to forget about the conditions.  Here is how this is done:

Integrate[1/((a-x)*(b-x)),{x,0,xt},GenerateConditions->False]

If you do not do that you get the following monster

ww=Integrate[1/((a-x)*(b-x)),{x,0,xt} ]

You do not know enough Mathematica programming to interpret this expression. You will learn more about the If function later and then you could use this expression.

Cell 9:   Numerical integration.

Sometimes Mathematica cannot perform an integral. In such cases it returns the expression unevaluated.

Integrate[Log[Sin[x]]/Sqrt[x],x]

It took me a while to find an integral the computer could not do. The one given above did it. When this happens you can perform the integral numerically.

NIntegrate[Log[Sin[x]]/Sqrt[x],{x,1,3}]

Numerical integration is performed only on a definite integral and the function to be integrated must contain no symbols other than the integration variable. If a symbol is present (a in the expression below) you get an error message and the computer returns the input unevaluated

Clear[a]
NIntegrate[a Log[Sin[x]]/Sqrt[x],{x,1,3}]

Sometimes the function NIntegrate provided by the machine fails.  In that case you will have to understand why the integral is difficult and use the methods describe in a book of numerical analysis to write a program that overcomes that difficulty.  Mathematica book offers a number of suggestions.

?Cell 10: How to solve equations analytically

Mathematica solves algebraic equation, especially if they are polynomials. The syntax is

	Solve[equation, unknownQuantity]
	
The unknown quantity is a symbol. Mathematica has a special way of writing equations as you can see below:

sol=Solve[a x^2+b x +c ==0,x]

Note the syntax defining the equation: the equality in an equation is == (two equal signs without a blank between them ) not = (the later sign is used for assignments only). The output is a list of replacement rules. In the case above there are two rules because the equation has two solutions. I assigned the result of the evaluation to the symbol sol.

The fact that the output is a list of lists can cause slight complications in the way we use it.

Let us explore how we can use the result. Let us say that we have a function f[x] and we want to evaluate it for the values of x given by the roots of the equation. You can do the following:

f[x_]:=x^2
fx=f[x]/.sol

I defined the function f[x] and then I used 

 	f[x]/.sol
 	
 to replace x in the function with the solutions of the equation. Since there are two solutions I get two results, which are returned as the elements of a list. I called the output list fx.

Whenever you can, it is better to use the syntax

f[x]/.Solve[a x^2 +b x +c ==0,x]

when you want to apply the function f to the roots of the equation a (x^2)+ b x +c =0. In this way you do not define the intermediate symbol sol. It is always best to have the smallest number of assignments and symbols.

Because the output

 	{{x->((-b-(Sqrt[(b^2)-4 a c]))/(2 a))},{x->((-b+(Sqrt[(b^2)-4 a c]))/(2 a))}}
 	
 is a list of two lists, each containing a replacement rule, the result  of 
 
 f[x]/.Solve[a x^2 +b x +c ==0,x]
 
 is as a list containing two values of f[x], with x replaced with one of the roots. Sometime you need to use the result in further calculations and having it in a list is inconvenient. Here are several ways of handling this

Method 1:  You can assign the two results to the symbols fx1 and fx2 as follows.

fx1=fx[[1]]

fx2=fx[[2]]

Here fx[[1]] means the first element of the list fx, etc.

Method 2: Use

Solve[a x^2 +b x +c ==0,x]

Solve[a x^2 +b x +c ==0,x]//Flatten

{f1,f2}=f[x]/.Solve[a x^2 +b x +c ==0,x]//Flatten

This automatically assigns the result to f1 and f2

f1

f2

The function Flatten removes the inner brackets from a list:

{{a,b,c},{A,B}}//Flatten

Mathematica also solves systems of equations. The syntax is 

	Solve[{eq1,eq2,...},{Unknown1, Unknown2,....}]
	
eq1, eq2, .... must be equations, which means that they have the form:  

	expression==anotherExpression

Clear[a,b,c,d,e,f,x,y]
sol=Solve[{a x+b y==c, d x +e y==f},{x,y}]//Flatten

If you want to use the value of x and y as arguments in the function g[x,y] use

g[x,y]/.(Solve[{a x+b y==c, d x +e y==f},{x,y}]//Flatten)

Another example:

x^2+y^2/.(Solve[{a x+b y==c, d x +e y==f},{x,y}]//Flatten)

This one-line program solves the two equations and than replaces x and y, in the expression (x^2) + (y^2), with the solutions

In most cases you do not want a general solution, but you want to solve the equations for specific numerical coefficients.

Clear[x,y,z];
eqs={3.2 x +2.1  y+3.7  z^2==2., 5.4  x+2.3  y^2 +3.7  z==2.1,2.1   x+ 7.  y +4.9  z==7.1};
sol=Solve[eqs,{x,y,z}]

There is more than one solution, because the equations are of second order. In this case you do not want to Flatten the output since different lists give you different solutions. You can pick up the first solution by using sol[[1]]

m[x,y,z]/.sol[[1]]

Or you can get a list of the results corresponding to all solutions:

m[x,y,z]/.sol

Cell 11: Numerical solutions of equations:  FindRoot and NSolve

Sometimes Solve is unable to provide a solution. In that case you can use a numerical procedure. This can be done only if all the coefficients and the exponents in the equation are numerical. No symbols, other than the unknown quantities, are allowed in the equation. The syntax is:

	FindRoot[equation, {x,(Subscript[x, 0])}]
	
Here x is the unknown and (Subscript[x, 0]) is a guess of the solution.  This does not have to be a good guess, but the  better the guess the more efficient the procedure for finding the solution.  Here is an example:

Clear[x]
eq=2.3 x^6+4.2 x^2 ==3.2;
sol=FindRoot[eq,{x,2.}]

I had no idea what the root was, so my guess was random. It did work.

Check that I did get a root.

2.3 x^6+4.2 x^2==3.2/.FindRoot[eq,{x,2.}]

The syntax   number1==number2 ( here == is two equal signs with no space between them) tests whether number1 is the same as number2. If this is the case the output is True. Otherwise it is False.

The equation

	2.3 x^6+4.2 x^2 ==3.2
	
is of sixth order and it should have six roots. One is x=0.791803.  To see what the others might be, make a plot of

	2.3 x^6+4.2 x^2-3.2
	
and see where it cuts the x axis (that is where the roots are)

Plot[2.3 x^6+4.2 x^2-3.2,{x,-1,1}];

There is a root at 0.79 and another one at -0.79 (as far as I can read it on the graph). I can find the negative one as follows

sol=FindRoot[eq,{x,-2.} ]

If I gave FindRoot a guess closer to the negative root, the computer finds the negative root

Here is a mystery for you: a sixth order equation must have six roots. We found two. Where have the other four roots disappeared? 
They must be complex numbers and FindRoot did not find them.

To find complex roots you have to give as starting point a complex number:

sol=FindRoot[eq,{x,-2.+I 3.} ]

I can check whether this satisfies the equation:

2.3 x^6+4.2 x^2==3.2/.FindRoot[eq,{x,-2.+I 3.} ]

Let us not hurry and conclude that the root is erroneous. The root is approximate and it may be that the evaluation of

	 2.3 (x^6) + 4.2 (x^2), 
when 
	x = 0.766879 + i 0.949515 
is close to 3.2 but not sufficiently close for the computer to declare them equal. Let's evaluate 2.3 (x^6)-4.2 (x^(2 ))ourselves and see if it is equal to 3.2:

2.3 x^6+4.2 x^2/.FindRoot[eq,{x,-2.+I 3.} ]

This is pretty close to zero 3.2, so the computer has found one complex solution.

Exercise:   Find the other complex roots of the equation studied above

FindRoot uses an iterative procedure and sometimes the solution requires a lot of iterations.  By default FindRoot has a limit on the number of iterations. We can find what this limit is by using:

Options[FindRoot]

Mathematica 5.1 allows 100 iterations.  If the initial guess is very bad we may need more. You can increase the number of iteration as shown below:

sol=FindRoot[eq,{x,2.}, MaxIterations->600]

Obviously in this case it was not necessary to have more iterations, since the solution for 600 iterations is the same as the solution for 100.

Polynomial equations are solved numerically by the function NSolve. This has the syntax

	NSolve[equation, unKnown]

Let us solve the equation eq define earlier. This is:

eq

sol=NSolve[eq,x]

Note that NSolve, unlike FindRoot, gives all six solutions as replacement rules. Inserting the roots in the expression 

	2.3 x^6+4.2 x^2-3.2
	
shows (see below) that all six roots given by NSolve provide a very good approximation to the six solutions. 
I check below that they are solutions:

2.3 x^6+4.2 x^2-3.2/.NSolve[eq,x]

If you got tired of seeing extremely small numbers in the result, you can remove them by using the function Chop

2.3 x^6+4.2 x^2-3.2/.NSolve[eq,x]//Chop

NSolve and FindRoot can be used to solve systems of equations. Here is how this is done:

Clear[x,y]
eqs={3 x^3+2 y^2==3,4 y^4-3 x==7}
sol1=FindRoot[eqs,{x,1},{y,2}]

I can check whether this is a solution by using:

{3 (x^3)+2 (y^2)-3,-3 x+4 (y^4)-7}/.FindRoot[eqs,{x,1},{y,2}]//Chop

Since both equations are polynomials NSolve finds all the roots:

sol2=NSolve[eqs,{x,y}]

Check that they are solutions:

{3 (x^3)+2 (y^2)-3,-3 x+4 (y^4)-7}/.sol2//Chop

Cell 12: A numerical method for finding the minimum of a function:

The function FindMinimum can be used for a numerical calculation of the minimum of a function.  I will use the function 

	3.2 x -x^3
	
I plot it to see where the minimum is:

Plot[3.2 x -x^3,{x,-10,10}]

If you click on the graph, hold Ctrl key down, and place the cursor on a point on the graph, you can read its coordinates at   on the left side of the lower bar in the file. I get x = -1.1 for the position of the minimum. This value is approximate since the reading of the graph is approximate

Let's first find the minimum by conventional method:
	solve    (df/dx)=0

sol=Solve[D[3.2 x -x^3,x]==0,x]

From the graph I know that x = -1.0328 corresponds to the minimum and x = 1.0328 corresponds to the maximum. To check algebraically that this is correct I calculate the second derivative at these two points

D[3.2 x -x^3,{x,2}]/.Solve[D[3.2 x -x^3,x]==0,x]

The first number in the output is the second derivative at the point x = -1.0328. Since this is positive, the point x = 1.0328 is a minimum. The second derivative at the other point is negative and that point corresponds to a maximum.

Start the search for the minimum, with FindMinimum, with the guess x = - 2.

FindMinimum[3.2 x -x^3,{x,-.2}]

The first number in the output is the minimum value of the function. The second term in the list is a list {x->number} where number is the value of x for which the function has a minimum. But what about the warning message?  This neurosis, a tendency to give an unreasonable number of warnings, is a new feature of Mathematica 5. The message tells us that the computer could not find the minimum with the accuracy required of it. Here is how you tell the computer that you would be happy with an accuracy at the fourth decimal point

FindMinimum[3.2 x -x^3,{x,-.2},AccuracyGoal->4]

Increasing the accuracy does not change the result, so we accept it.

FindMinimum[3.2 x -x^3,{x,-.2},AccuracyGoal->5]

FindMinimum[3.2 x -x^3,{x,-.2},AccuracyGoal->6]

If I started with another guess I could have had trouble.

FindMinimum[2.3 x-x^3,{x,2}]

Something is awfully wrong: overflow means that the numbers are becoming larger than what the computer can handle. Look at the graph to see why. I started at the right of the right-most peak, and the program follows the function down hill and finds no minimum in maximum number of iterations. This is as it should be, since there is no minimum in the region of the search. The search keeps going down hill making x larger and larger and it stops when the maximum number of iterations (which is 100) was reached.

If I start with x=0, the minimum is found:

FindMinimum[3.2  x-x^3,{x,0.}]

FindMinimum will work for functions of many variables. Let us try first an easy function, which has a minimum at x=0, y=0 , z=0

Clear[x,y,z]
FindMinimum[x^2 +y^2+z^2,{x,0.1},{y,.2},{z,1.}]

Here is another function of three variables

FindMinimum[2.1 x^3+.1  (x-y) y +y^2+z^2+ 2.,{x,0.1},{y,.2},{z,1.}]

If a function has several minima, FindMinimum tends to find the minimum located downhill from the starting point (remember that iterative methods for finding minima move down hill from the starting point). Here is an example:

Plot[(x-2)^3 Sin[3 x],{x,-6,6}]

Note that which minimum is found depends on the starting point:

FindMinimum[(x-2)^3 Sin[3 x],{x,-0.5}]

FindMinimum[(x-2)^3 Sin[3 x],{x,-1.}]

FindMinimum[(x-2)^3 Sin[3 x],{x,-0.65}]

FindMinimum[(x-2)^3 Sin[3 x],{x,-4.}]

Cell 23: Solving differential equations

I will solve first a simple differential equation for the harmonic oscillator:

	(((d^2)x(t))/(dt^2)) =- (ω^2)x(t)			(1)
	
with the initial conditions 

	x(0)=x0				(2)
	
and
	((dx(t))/dt) (Subscript[|, t=0]) =v0				(3)

I will use the function DSolve. The syntax is 

	DSolve[equation, x[t],t]
	
The first argument is the differential equation, the second argument indicates the function I solve for, the third arguments indicates that t is the argument in the function I solve for. For my example the equation is:

	 D[x[t],{t,2}] == -ω^2 x[t];
	 
an equivalent (and simpler) way of writing the equation is:

	x''[t]== -ω x[t]
	
Since I did not impose the initial conditions, the solution will contain two arbitrary constants, for which Mathematica reserved the symbols C[1] and C[2].

The equation to be solved is:

Clear[x]
eq=x'' [t] ==-ω^2 x[t]

Note that when I wrote the differential equation I have used x''[t] for the second derivative. The first derivative is x'[t], the third is x'''[t].  I could have used D[x[t],{t,2}], but this involves more typing and it is rarely used. Let us see how this notation works

Clear[f]
f[q_]:=q^2
f'[q]

f''[r]

f'''[q]

One more example:

Clear[f]
f[v_]:=v^2 Cos[v]
f'[c]

f''[d]

You can take more examples if you want to convince yourself that the syntax f '[x] generates the first derivative of f with x, 
f ''[x] gives the second derivative, etc

Let us get back to the differential equation and solve it:

eq

Clear[x]
sol=DSolve[ eq,x[t],t]

The solution is given as a replacement rule. Since there is  only one solution I do not need a list of a list with one element. I flatten the list sol and assign the result to sol. This redefines the meaning of sol

sol=sol//Flatten

Create a new function f which is the solution of the differential equation:

Clear[f]
f[z_]:=(x[t]/.sol)/.t->z

Clear[f]
f[t_]:=x[t]/.sol

f[z]

Test that the function f works as intended

f[x]

Above I use the replacement rule in  x[t]/.sol to replace x[t]  with the x[t] given in the rule sol (which is the solution of the equation).  Then I use /.t->z to replace the argument t in the expression for x[t] with the argument z ( z is the argument in the definition of the function f[z_]:= ) .  This may seem strange and I'll explain later why you need to do it in this way.  As you have seen in the cell above, f defined in this way works

We have now created a function f[r] which is the solution of the differential equation 
	The line 
	
		f[z_]:=(x[t]/.sol)/.t->z  
		
is bad programing. The definition of the function assumes that the symbol sol will never be reassigned another value at any later time in the program. It is unsafe to assume this. If another value is assigned to sol, then the function f[z] will use it because the delayed assignment in  f[z_]:=(x[t]/.sol)/.t->z  means that the function is evaluated anew every time the symbol f[argument] is encountered in the program. Every time sol gets a new assignment, the function will use the new meaning of sol. It will no longer be the function you intended. One can avoid this problem by using the syntax shown below:

Clear[f]
f[r_]:=x[t]/.(DSolve[x'' [t] ==-ω^2 x[t],x[t],t]//Flatten)/.t->r

f[z]

There is now no sol in the syntax and we no longer can ruin the calculation of f[x] by making another assignment to sol.  However, this syntax is also troublesome, because using the delayed assignment := the differential equation is solved anew every time f[x] is used. If the function f[x] is used many times it is best to use the immediate assignment =  (see below)

Clear[g]
g[r_]=x[t]/.(DSolve[x'' [t] ==-ω^2 x[t],x[t],t]//Flatten)/.t->r

I have not used a delayed assignment :=, but the assignment =   This means that the expression at the right is evaluated only once, when it is defined. Had I used := the function would have been evaluated every time the expression f[r] was encountered in the program. That means that every time f[r] is encountered the computer solves the differential equation again! You can test how much time is wasted in this way

We can check the difference between f and g:

?f

?g

f stores the program that calculates the solution of the differential equation, while g stores the solution. You see now, in action, the difference between := and =

Determine the integration constants: The solution g[t] depends on the integration constants C[1] and C[2]. I determine their value by forcing the solution to satisfy the initial conditions 
	
	g[t=0]=x0 
and
	g'[t=0]=p0

Here are the equations corresponding to the two initial conditions:

eqs={g[0]==x0,(g'[t]/.t->0) ==p0}

Note that

g'[0]

Solve these equations for C[1] and C[2]

Solve[eqs,{C[1],C[2]}]//Flatten

I use this  result to replace C[1] and C[2] in the solution g[t] of the differential equation.

solution[x_]=(g[t]/.(Solve[eqs,{C[1],C[2]}]//Flatten))/.t->x

Note the use of = (not :=) and that I do not give the solution of the equation a name (such as sol in the earlier cells). This function works

g is redefined by the above statement. Let's check that it is what we want.

?solution

Clear[z]
solution[z]

This is the solution of the differential equation (1) with the boundary conditions (2) and (3).

Exercises:  Solve the following differential equations:

y'[x]=(2.+x y[x]) with initial condition  y[0]=1.

y''[t]+.2 x y'[t]+3. y[t]==0  with y[0]=1. and y'[0]=3.1

Write a few equations yourself and see if you can solve them. 
In all cases create a function that corresponds to the solutions.