Using Cubic Splines in Practice

Posted on 06/05/2010 by dougaj4

I recently had an interesting query from Georg,  a geography student in Germany regarding the use of cubic splines to fit a smooth curve to experimental data.  The graph he sent me is shown below:

The cubic spline has produced a nice smooth curve, but the problem is that it is supposed to be a cumulative frequency curve, and cumulative frequency curves are not supposed to slope downwards, let along pass below the zero line.  Unfortunately cubic splines don’t know that, and will just adopt the shape dictated bt the given points.

One solution is to insert additional points in the areas with problems as illustrated below:

Original data

Insert two additional points

 Which gives us the result we want!

The new all ascending curve

Posted in Charts, Excel, Maths, Newton, UDFs | Tagged , | 9 Comments

Numerical solutions with CSpline

Posted on 05/05/2010 by dougaj4

The CSpline function presented in a previous post  fits a series of cubic polynomial curves to a specified series of points, returning the Y values for listed intermediate X values.  I was recently asked if this could be reversed to find the X value for a specified Y value.  Since CSpline can also return the coefficients defining each of the cubic curves it turns out that this is fairly simple.  The procedure is:

  1. Use CSpline to find the coefficients defining each segment of the spline, that is find the value of a, b, c, and d in Y = aX^3 + bX^2 + cX + d.
  2. For each Y value:
    1. For each segment
      1. Solve the applicable cubic equation for the specified Y value
      2. Check if any solutions lie within the X range of the current segment
      3. If a solution is found add, it to the results array and go to the next Y
      4. If all solutions are outside the segment X range go to the next segment
    2. If no solution is found for any segment return a message to the results array and go to the next Y
  3. Assign the results array to the function return value.

The new version of CSpline.xls, including the new function, SolveSplineA, can be downloaded from: CSpline2.zip

Some screenshots of examples are shown below:

Solvespline input and output (click for full size view)

Example 1

Examples 2 and 3

Posted in Excel, Maths, Newton, UDFs, VBA | Tagged , , , , | Leave a comment

More about not being different

Posted on 29/04/2010 by dougaj4

Last year I posted a link to a puzzle at Tanya Khovanova’s Math Blog.  I just came across a follow up to that puzzle:

The Odder One Out

where she discusses the responses. 

Which by some devious path reminded me of my very favourite bit of Monty Python’s Life of Brian:

Posted in Bach | Tagged , | Leave a comment

Automating chart scale limits

Posted on 28/04/2010 by dougaj4

Edit 22 Mar 2014:  Also see https://newtonexcelbach.wordpress.com/2012/11/25/automating-chart-scale-limits-update/
for improved version with download link, and example of plotting a chart from a formula entered as text.

One of the more annoying things about Excel charts is that if you want to over-ride the automatic scale limits there is no built-in way to link the limits to a spreadsheet cell, so you have to go into the dialog box and change the numbers manually every time you want to change the scale.

Tushar Mehta recently posted a solution at Daily Dose of Excel, with an add-in providing this functionality, as well as other format chart adjustments.

In the following discussion John Walkenbach linked to an earlier solution that those who like to keep things simple may prefer.  This is a user defined function (UDF), that returns no data, but modifies the y-axis limits of any named chart.  It should be noted that this solution only works in Excel 2007 and later, is undocumented, and is not supposed to work at all, so use with caution.  I have taken the liberty of modifying John’s code so that the limits of both the X and Y axes can be linked to cell values, the axes limits can be re-set to automatic, and the status of each axis is returned by the function.  The revised code and a screen shot are given below.

Function ChangeChartAxisScale(CName As String, Optional Xlower As Double = 0, Optional Xupper As Double = 0, _
Optional Ylower As Double = 0, Optional YUpper As Double = 0) As Variant
Dim Rtn As String
'   Excel 2007 only

With ActiveSheet.Shapes(CName).Chart.Axes(1)
If Xlower = 0 Then
.MinimumScaleIsAuto = True
Rtn = "Xmin = auto"
Else:
.MinimumScale = Xlower
Rtn = "Xmin = " & Xlower
End If

If Xupper = 0 Or (Xupper < Xlower) Then
.MaximumScaleIsAuto = True
Rtn = Rtn & "; Xmax = auto"
Else
.MaximumScale = Xupper
Rtn = Rtn & "; Xmax = " & Xupper
End If
End With

With ActiveSheet.Shapes(CName).Chart.Axes(2)
If Ylower = 0 Then
.MinimumScaleIsAuto = True
Rtn = Rtn & "; Ymin = auto"
Else
.MinimumScale = Ylower
Rtn = Rtn & "; Ymin = " & Ylower
End If

If YUpper = 0 Or (Xupper < Xlower) Then
.MaximumScaleIsAuto = True
Rtn = Rtn & "; Ymax = auto"
Else
.MaximumScale = YUpper
Rtn = Rtn & "; Ymax = " & YUpper
End If
End With

ChangeChartAxisScale = Rtn
End Function

ChangeChartAxisScale Function, Click for full size view

Posted in Charts, Excel, UDFs, VBA | Tagged , , , , | 8 Comments

Analysis of RC circular sections – theory and code

Posted on 25/04/2010 by dougaj4

The basic tasks involved in analysing the stresses and strains in a reinforced concrete section under combined bending and axial load are:

Assuming a horizontal beam with the top of the beam in compression and the bottom in tension

  1. Find the depth of the Neutral Axis below the top face, that is the depth of the line along which the strain is zero
  2. Find the stress in the concrete at the top face
  3. Find stresses at other levels by simple proportion
  4. Find forces and bending moments
  5. Confirm that the calculated total force and moment in the section are equal to the applied load

The depth of the Neutral Axis is controlled by the eccentricity of the reaction force, which must be equal to the eccentricity of the applied load:

Mr/Pr = M/P = E

where the bending moment is taken about any convenient axis, which for the purposes of this analysis will be taken as the top fibre of the section.

It can be shown that:

Pr= KQx
and
Mr = KIx

where Qx and Ix are the first and second moments of area about the Neutral Axis, at depth X below the top (most compressive) fibre of the section, and K = Top fibre stress / X

Hence  KIx / KQx = Ix / Qx = E + X

Polynomial equations to determine the value of X are derived for any section made up of a series of trapezoidal layers in a previous article, but if we wish to apply this equation to a circular section it is more convenient to solve by iteration.  The actual equation used in the analysis presented here is:

Qx(E + X) / Ix – 1 = 0

which behaves nicely as Q approaches zero.

The procedure to find X for a circular section is then:

  1. Assuming that the reinforcement is distributed evenly around the circle, find the depth of each reinforcement layer from the top face.
  2. Find the transformed area of a reinforcing bar in tension and compression.  That is the bar area x Steel Elastic Modulus / Concrete Elastic Modulus in tension, (since the bar displaces an equal area of concrete) bar area x (Steel Elastic Modulus / Concrete Elastic Modulus – 1) in compression.
  3. Find the transformed area properties (area, level of centroid, first moment of area about layer depth,   of area about centroid depth) for the neutral axis at each steel layer, adjusting the transformed area for the numbers of bars in compression and tension.
  4. Find the composite (concrete + steel) Ix and Qx for the Neutral Axis at each reinforcement layer, and hence find which layers it lies between.
  5. Use Brent’s Method to find the Neutral Axis position that satisfies the governing equation given above.
  6. Find the composite area properties for the Neutral Axis at this layer.
  7. Hence find the stresses, strains, forces and bending moments in the section.
  8. Check that the calculated total axial force and bending moment are equal to the applied loads.

Excerpts of the VBA code are shown below.  Full open source code is included in the download file: RC design functions5.zip

Composite area properties about an axis at depth X below the top fibre:

Function CircleQxoI(SectData As Variant, X As Double) As Variant
Dim R As Double, H As Double, Theta As Double, A As Double, Xcr As Double
Dim Icc As Double, Iccent As Double, Icb As Double, Qcb As Double
Dim Ast As Double, Istc As Double, Xsc As Double, Ecc As Double
Dim Qsb As Double, Isb As Double, outA(1 To 8, 1 To 1), CallFunc As Long

  If TypeName(SectData) = "Range" Then SectData = SectData.Value2
  R = SectData(1)
  Ast = SectData(2)
  Xsc = SectData(3)
  Istc = SectData(4)
  Ecc = SectData(5)
  CallFunc = SectData(6)

    ' Concrete properties
    If X > 0 Then
        If X < 2 * R Then
            H = (R ^ 2 - (R - X) ^ 2) ^ 0.5
            Theta = WorksheetFunction.Atan2((R - X), H)
            A = R ^ 2 * (Theta - (Sin(2 * Theta) / 2))
  Xcr = (2 * R / 3) * ((Sin(Theta) ^ 3) / (Theta - Sin(Theta) * Cos(Theta)))
  Qcb = A * (X - (R - Xcr))
  Icc = A * R ^ 2 / 4 * (1 + (2 * (Sin(Theta)) ^ 3 * Cos(Theta)) / ((Theta - Sin(Theta) * Cos(Theta))))
  Iccent = Icc - A * Xcr ^ 2
        Else
            A = Pi * R ^ 2
  Qcb = A * (X - R)
  Iccent = A * R ^ 2 / 4
        End If
  Icb = Iccent + A * (Xcr - (R - X)) ^ 2

    Else
  Qcb = 0
  Icb = 0
    End If
    ' Reinforcement properties
  Qsb = Ast * (X - (R - Xsc))
  Isb = Istc + Ast * ((R - Xsc) - X) ^ 2

  If CallFunc = 1 Then
  CircleQxoI = ((Qcb + Qsb) * (Ecc + X) / (Icb + Isb)) - 1

    Else
        outA(1, 1) = X
  outA(2, 1) = Icb
  outA(3, 1) = Isb
        outA(4, 1) = A
        outA(5, 1) = Ast
        outA(6, 1) = A + Ast
  outA(7, 1) = Qcb
  outA(8, 1) = Qsb
  CircleQxoI = outA
    End If

End Function

Use Brent’s Method to find the depth of the Neutral Axis within a specified layer

Function QuadBrent(Func As String, ByVal Ak As Double, ByVal Bk As Double, Coefficients As Variant, _
  Optional YTol As Double = 0.00000000000001, Optional MaxIts As Long = gMaxIts, _
  Optional Subr As Long = 1, Optional XTol As Double = 0.0000000001) As Variant
  Dim Ck As Double, Sk As Double, Dk As Double, XRange As Double
  Dim Yak As Double, Ybk As Double, Yck As Double, Ysk As Double, Ydk As Double
  Dim i As Long, LastStepB As Boolean, Temp As Double, ResA(1 To 3, 1 To 1) As Double, SubName As String

  If TypeName(Coefficients) = "Range" Then Coefficients = Coefficients.Value2

  Yak = Application.Run(Func, Coefficients, Ak)
  Ybk = Application.Run(Func, Coefficients, Bk)

  If (Yak < 0 And Ybk < 0) Or (Yak > 0 And Ybk > 0) Then
  QuadBrent = "Yak and Ybk must have different signs"
        Exit Function
    End If

  If Abs(Yak) < Abs(Ybk) Then
        Temp = Ak
        Ak = Bk
        Bk = Temp
        Temp = Yak
  Yak = Ybk
  Ybk = Temp
    End If

    Ck = Ak
  Yck = Yak
  LastStepB = True

  For i = 1 To MaxIts

  XRange = (Bk - Ak)

  If (Abs(Ybk) < YTol) Then Exit For

  If Yak <> Yck And Ybk <> Yck Then

  ' SolveInvQuad uses inverse quadratic interpolation, SolveQuad uses Mullers method of quadratic interpolation
  If Subr = 1 Then SubName = "solveInvQuad" Else SubName = "solveQuad"
  Sk = Application.Run(SubName, Ak, Yak, Bk, Ybk, Ck, Yck) '

        Else
  Sk = Bk - (Ybk) * (Bk - Ak) / (Ybk - Yak)
        End If

        If ((Sk - (3 * Ak + Bk) / 4) * (Sk - Bk)) > 0 Then
            Sk = (Ak + Bk) / 2
  LastStepB = True
  ElseIf LastStepB = True Then
  If Abs(Sk - Bk) >= Abs(Bk - Ck) / 2 Or (Abs(Bk - Ck) < Abs(XTol)) Then
                Sk = (Ak + Bk) / 2
  LastStepB = True
            Else
  LastStepB = False
            End If
  ElseIf Abs(Sk - Bk) >= Abs(Ck - Dk) / 2 Or Abs(Ck - Dk) < Abs(XTol) Then
            Sk = (Ak + Bk) / 2
  LastStepB = True
        Else
  LastStepB = False
        End If

  ' Use Application.run to evaluate Func for the value Sk
  Ysk = Application.Run(Func, Coefficients, Sk)
        '----------------------------------------
        Dk = Ck
  Ydk = Yck
        Ck = Bk
  Yck = Ybk

  If Yak * Ysk < 0 Then
            Bk = Sk
  Ybk = Ysk
        Else
            Ak = Sk
  Yak = Ysk
        End If
  If Abs(Yak) < Abs(Ybk) Then
            Temp = Ak
            Ak = Bk
            Bk = Temp
            Temp = Yak
  Yak = Ybk
  Ybk = Temp
        End If

  Next i

  If Abs(Ybk) < Abs(Ysk) Then
  ResA(1, 1) = Bk
  ResA(3, 1) = Ybk
    Else:
  ResA(1, 1) = Sk
  ResA(3, 1) = Ysk
    End If
  ResA(2, 1) = i - 1

  QuadBrent = ResA
End Function
Posted in Beam Bending, Concrete, Excel, Newton, UDFs, VBA | Tagged , , , , , , | Leave a comment