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