MathGroup Archive 1995

[Date Index] [Thread Index] [Author Index]

Search the Archive

mathematica syntax, matrix notation

  • To: mathgroup at smc.vnet.net
  • Subject: [mg2127] mathematica syntax, matrix notation
  • From: itcook at peseta.ucdavis.edu (Poor Little Rich )
  • Date: Wed, 4 Oct 1995 01:55:24 -0400
  • Organization: University of California, Davis

Hello, I am a poor foolish undergrad in chemical engineering,and am 
working on a project which requires that I  understand the 
syntax of the following expression.  
The context is that I am trying to solve a system of ode's of the 2
point 
BV type.  That's probably not important to know.  What I would like to 
know is, what significance do the bracketed and double-bracketed 
expressions have in the following mathematica command:

yend[#189#_]:=(soln=NDSolve[system[#189#],{y1,y2},
                                                {x,-1,1}])[[1,1,2]][1]

It has been suggested to me (and I am paraphrasing) that the double
brackets indicate that the variable "soln" will be of a nested index
function form and that the double brackets reference a part of this
function, e.g:

soln={{{a,b},{a2,b2}},{{a3,b3},{a4,b4}}}

and thus passes a value of 1 to the function I have called "b"
here,then passes b[1] to yend[#189#_], as in:
yend[#189#_]:=(soln)[b[1]].
  
Does that seem correct?

A related question is that when I enter the following command:

system[#189#_]={y1'[x]==y2[x], y2'[x]==-Sin[y1[x]]+Cos[5 x],
                        y1[-1]==0,y2[-1]==#189#}

mathematica echoes back the same thing, just as you would expect,
except 
that it inexplicably adds a #1  right before the last #189#, like this:

system[#189#_]={y1'[x]==y2[x], y2'[x]==-Sin[y1[x]]+Cos[5 x],
                        y1[-1]==0,y2[-1]==#1 #189#}

Does anyone know why this is?  Thank you for your help in this matter. 


If it helps, I am enclosing the source of these statements, a
discussion 
on using mathematica to solve 2-point BVODE's:



Boundary Value Problems
Introduction
In this tutorial we will discuss how to use Mathematica to solve n-th
order
boundary value problems (BVP). Since we can always write a n-th order
BVP as
a system of n  first-order equations of the form y'(x)=F(y,x,u), where
y={y1,y2,..yn}, and F={F1,F2, F3,..Fn} and u={u1,u2,...} , our
discussion
will assume that the BVP to be solved can be written in this manner.

Mathematica has a numerical differential equation solver called
NDSolve[]
which solves initial value problems. Thus in order to use this routine
for
BVPs, we must implement some iteration procedure that is based on a
guess for
the unknown (missing) initial condition. One way to accomplish this is
to use
a shooting method. Suppose we have a 2nd-order BVP for y(x) subject to
the
BCs y1(a)=b1, y2(b)=b2. To integrate this system from x=a to x=b, we
need to
guess the value of y2(a). Since our guess is unlikely to be correct,
the
numerical solution y(x) will not satisfy the BC y2(b)=b2. Thus we must
update
our guess for y2(a), and integrate once again.

Since the numerical solution y(x) depends parametrically on the initial
guess
y2(a)=#189#, we have an additional equation for #189#, viz., y2(x;
#189#)=b2.
Note that we have expressed the dependence on #189# as y(x; #189#). To
find
#189# we augment the system of ODEs with the additional equation, and
solve
iteratively for #189#. We now demonstate this procedure in the next
section
 
Example 1:
We consider the following BVP defined by the differential equation
y''(x)+sin(y(x))=cos(5x)
with boundary conditions
y(-1)=0,  y(1)=0
This equation can be expressed as a system of two first order ODEs by
making
the substitutions:
y1(x)=y(x)
y1'(x)=y2(x)
y2'(x)=-sin(y(x)+cos(5x)
with boundary conditions
y1(-1)=0,   y2(-1)=#189#
Next we assign the system of equations to a variable system:
system[#189#_]={y1'[x]==y2[x], y2'[x]==-Sin[y1[x]]+Cos[5 x],
                        y1[-1]==0,y2[-1]==#189#}
Note that the variable system is a function of the unknown initial
condition
#189#. Next, we define a function yend[#189#_] which is equal to the
solution
of the ODEs at the end point x=1.Refer to the tutorial on the use of
NDSolve
for an explanation of the format used in this statement 
yend[#189#_]:=(soln=NDSolve[system[#189#],{y1,y2},
                                                {x,-1,1}])[[1,1,2]][1]
We use the Plot[] function to show how the end point yend[#189#_]
varies with
the initial condition
Plot[yend[#189#],{#189#,-2,2}]
The value of #189# that we are seeking is the value that satisfies the
equation yend[#189#]==0. We make use of the FindRoot[] function to
determine
#189# such that yend[#189#]= y1(1)=0:
FindRoot[yend[#189#]==0,{#189#,-2,2}][[1,2]]
Finally, we use the Plot[] routine to plot out the solution y1[x].
Plot[Evaluate[y1[x]/.soln],{x,-1,1}]
Example 2:
In this example we show how to use more advanced programming in
Mathematica    
to combine all the steps in Example 1 so as to have an efficient
algorithm.
The first several ines of the code are similar except that we make use
of
"pure functions" defined using the symbols # and &:
system[#189#_]={y1'[x]==y2[x], y2'[x]==-Sin[y1[x]]+Cos[5 x],
                                                                       
                                y1[-1]==0,y2[-1]==#189#}
yend:=(soln=NDSolve[system[#],{y1,y2},
                                                       
{x,-1,1}])[[1,1,2]][1]&
FindRoot[yend[#189#]==0,{#189#,-2,2}][[1,2]]
Plot[Evaluate[y1[x]/.soln],{x,-1,1}]
Eample 3:      
In this example we show how the ideas we have developed in the previous
sections can be applied to more difficult BVPs. The first BVP we solve
is the
3rd-order Blasius equation that describes the velocity profile f(n) in
a
laminar boundary layer over a flat plate:
ff''+2f'''=0,
f(0)=0,  f'(0)=0,  f'(#176#)=1
 As before, we reduce the 3rd-order equation to a system of 3
first-order
equation with the following substitutions:
 f'(n)=v(n),   
         v'(n)=w(n)
         w'(n)=-f(n)*w(n)/2
         with initial conditions
 f(0)=0,   v(0)=0,  w(0)=#189#
 and the auxillary equation to find #189#
         v(#176#; #189#)=1     
 In the calculations we truncate the infinite domain and apply the BC
at
infinity at a finite location given by n=6

system[#189#_]={2 w'[n]+f[n] w[n]==0,w[n]==v'[n],v[n]==f'[n],
                        w[0]==#189#,v[0]==0,f[0]==0}
gun:=(soln=NDSolve[system[#],{w,v,f},
                                        {n,0,10}])[[1,2,2]][10.]&
FindRoot[gun[#189#]==1,{#189#,1,0.5}]
Plot[Evaluate[{f[n],v[n],w[n]}/.soln],{n,0,10}]
In our next BVP we consider the differential equation for stagnation in
plane
flow (Hiemenz flow). The similarity form for the velocity field is
given by
the following third-order equation:
f'''+ff''-f'2+1=0
with boundary conditions
f(0)=0,  f'(0)=0,  f'(#176#)=1
We reduce this to a system of 3 first-order equations with the
following
substitutions:
 f'(n)=v(n),   
         v'(n)=w(n)
         w'(n)=v(n)2-1-f(n)w(n)         
with initial conditions
         f(0)=0,   v(0)=0,  w(0)=#189#  
 and the auxillary equation to find #189#
 v(#176#; #189#)=1
system[#189#_]={w'[n]==v[n]^2-1-f[n] w[n],f'[n]==v[n],
                        v'[n]==w[n],f[0]==0,v[0]==0,w[0]==#189#}
gun:=(sol=NDSolve[system[#],{f,v,w},{n,0,5}])[[1,2,2]][5]&
FindRoot[gun[#189#]==1,{#189#,-2,2}]
Plot[Evaluate[{f[n],v[n],w[n]}/.sol],{n,0,5}]  
                   2
{w'[n] == -1 + v[n]  - f[n] w[n], f'[n] == v[n], 
 
  v'[n] == w[n], f[0] == 0, v[0] == 0, w[0] == #189#}
{#189# -> 1.23259}

-Graphics-
^*)


---Sincerely, Rich Cook |"The volume of a body of rotation is equal to
the
  rdcook at ucdavis.edu    |generating area times the distance traveled by
the 
  1-916-758-6897        |centroid of the area while the body is being
  Metaphor Alert        |generated."         -Beer, et al.


-Sincerely, Rich Cook |"There are trivial truths & there are great 
rdcook at ucdavis.edu    |truths. The opposite of a trivial truth is 
1-916-758-6897        |plainly false. The opposite of a great truth    
                Metaphor Alert        | is also true."           -Neils
Bohr


  • Prev by Date: Re: Converting table to a number?
  • Next by Date: A simple swap function
  • Previous by thread: RE:MmaToHTML_for_Macs
  • Next by thread: Q: Graphical Features.