Keeping track of the Chain Rule is getting annoying

Lets do the computation on a general interval and that way we can skip that pain.

Interval [a,b] Computation

Note the Determinant of the matrix is (((a-b))^4) which is non-zero for interesting intervals.

Clear[basis,a,b,x]
cubics[x_]={1,x,(x^2),(x^3)};
abIntMatrix={
cubics[a],
cubics[b],
cubics'[a],
cubics'[b]
}
basis[{a_,b_}][x_]=FullSimplify[Transpose[Inverse[abIntMatrix]]].cubics[x];
TableForm[basis[{a,b}][x]]

({{(((3 a-b) (b^2))/(((a-b))^3))-((6 a b x)/(((a-b))^3))+((3 (a+b) (x^2))/(((a-b))^3))+((2 (x^3))/(((-a+b))^3))}, {(((a^2) (a-3 b))/(((a-b))^3))+((6 a b x)/(((a-b))^3))-((3 (a+b) (x^2))/(((a-b))^3))+((2 (x^3))/(((a-b))^3))}, {-((a (b^2))/(((a-b))^2))+((b (2 a+b) x)/(((a-b))^2))-(((a+2 b) (x^2))/(((a-b))^2))+((x^3)/(((a-b))^2))}, {-(((a^2) b)/(((a-b))^2))+((a (a+2 b) x)/(((a-b))^2))-(((2 a+b) (x^2))/(((a-b))^2))+((x^3)/(((a-b))^2))}})!margaret

Lets check this works.

{a,b}={0.1,0.2};
Plot[
basis[{a,b}][x],
{x,a,b}]

So a general "slab" function with specified value and slope at x=a and which is zero and flat at x=b is

{a,b}={0.1,0.3};
{v,s}={1.2,3.4};
slab[x_]=v basis[{a,b}][x][[1]] + s basis[{a,b}][x][[3]];
Plot[slab[x],{x,a,b},
Epilog->{Red,Map[ Point, {{a,v},{a,s},{b,0}}]}
]

To make a bump function we need to put two "matched" slabs together.

Clear[slab,x]
{xp,xc,xn}={0.1,0.3,0.5};
{v,s}={1.2,3.4};
slab[x_]=Piecewise[{
{v basis[{xc,xp}][x][[1]] + s basis[{xc,xp}][x][[3]],xp≤x≤xc},
{v basis[{xc,xn}][x][[1]] + s basis[{xc,xn}][x][[3]],xc≤x≤xn},
{0, True}
}
];
Plot[slab[x],{x,xp,xn},
Epilog->{Red,Map[ Point, {{xc,v},{xc,s},{xp,0},{xn,0}}]}
]

Putting it together. Interpolation

Tidying up the computation up above and calling the pieces we need vϕ and sϕ (for value ans slope ϕs).  The fisrt thing to do is make sure that we have everything put together correctly.

Clear[vϕ,sϕ,a,b,x]
vϕ[{a_,b_}][x_]:= (((3 a-b-2 x)(((b-x))^2))/( (((a-b))^3))) 
sϕ[{a_,b_}][x_]:= (((x-a) (((b-x))^2))/(((a-b))^2))
{a,b}=RandomReal[{0,1},2];
{va,sa}={1.2,3.4};
cubic[x_]=va vϕ[{a,b}][x]+sa  sϕ[{a,b}][x];
Plot[ 
cubic[x],
{x, Min[a,b],Max[a,b]},
Epilog->{
{Red, Line[{{a,0},{a,10}}]},
{Green, Line[{{b,0},{b,10}}]}
}
]

Now lets see what our approximations look like on a 1D mesh. Here is the interval/element based computation.

GridPts={0.0,0.1,0.2, 0.3, 0.5,1.0};
ValuesAndSlopes={{0.,0.},{0.4,1.4},{0.8,1.8},{1.2,2.2},{-3.4,2.4},{0.,0.}};

n=Length[Values];
pwCubic[x_]= Piecewise[
Table[
{a,b}=GridPts[[{i,i+1}]];
{{va,sa},{vb,sb}}=ValuesAndSlopes[[{i,i+1}]];
{va vϕ[{a,b}][x] + vb vϕ[{b,a}][x]+sa sϕ[{a,b}][x] + sb sϕ[{b,a}][x], a<x<b},
{i, 1, n-1}
]
];
Plot[pwCubic[x],{x,0,1},
Epilog->{Red, Map[Point,Transpose[{GridPts,ValuesAndSlopes[[All,1]]}]]}
]
Plot[pwCubic'[x],{x,0,1},
Epilog->{Red, Map[Point,Transpose[{GridPts,ValuesAndSlopes[[All,2]]}]]}
]

Here is the grid/node based computation

GridPts={0.0,0.1,0.2, 0.3, 0.5,1.0};
ValuesAndSlopes={{0.,0.},{0.4,1.4},{0.8,1.8},{1.2,2.2},{-3.4,2.4},{0.,0.}};

n=Length[Values];
pwCubicNode[x_]= Sum[
{xp,xc,xn}=GridPts[[{i-1,i,i+1}]];
{v,s}=ValuesAndSlopes[[i]];
Piecewise[{
{v vϕ[{xc,xp}][x] +s sϕ[{xc,xp}][x] , xp≤x≤xc},
{v vϕ[{xc,xn}][x] +s sϕ[{xc,xn}][x] , xc≤x≤xn}
}],
{i, 2, n-1}
];
Plot[pwCubicNode[x],{x,0,1},
Epilog->{Red, Map[Point,Transpose[{GridPts,ValuesAndSlopes[[All,1]]}]]}
]
Plot[pwCubicNode'[x],{x,0,1},
Epilog->{Red, Map[Point,Transpose[{GridPts,ValuesAndSlopes[[All,2]]}]]}
]

Putting it together. FEA

Now we have our functions correct we can go and build our stiffness matrix for the beam equation
	(Subscript[u, xxxx])=f
which has weak (or variational form
	(u'',v'')=(f,v)  for all v in a suitable space.
We substitute a sum of our shape functions into the weak form and get an equation multiply by each of our shape functions. we can do this in either a node based or an element based way.

Just to be completely clear at each grid point (Subscript[x, i]) (or node in the books other terminology) there are two shape functions:  a value shape function which gives a value 1 and slope zero at (Subscript[x, i]) and zeros at all the other nodes; and a slope shape function that gives value 0 and slope 1 at (Subscript[x, i]) and zeros at all the other nodes. Just like in the 1D piecewise linear case each of these shape functions is made up of two pieces one on each side i.e. [(Subscript[x, i-1]),(Subscript[x, i])]  and [(Subscript[x, i]),(Subscript[x, i+1])].  since we built our shape functions generally we can use the same recipe on each side.

Node based computation

As before the node based computation is best suited for hand computation

Computation

Clear[a,b,x]
Piecewise[{
{Simplify[vϕ[{(Subscript[x, i]),(Subscript[x, i-1])}]''[x]], (Subscript[x, i-1])<x<(Subscript[x, i])},
{Simplify[vϕ[{(Subscript[x, i]),(Subscript[x, i+1])}]''[x]], (Subscript[x, i])<x<(Subscript[x, i+1])}
}]//TraditionalForm

Clear[a,b,x]
Piecewise[{
{Simplify[sϕ[{(Subscript[x, i]),(Subscript[x, i-1])}]''[x]], (Subscript[x, i-1])<x<(Subscript[x, i])},
{Simplify[sϕ[{(Subscript[x, i]),(Subscript[x, i+1])}]''[x]], (Subscript[x, i])<x<(Subscript[x, i+1])}
}]//TraditionalForm

The value shape function on the intervals [(Subscript[x, i-1]),(Subscript[x, i])] and [(Subscript[x, i]),(Subscript[x, i+1])] is 
 	 (Subscript[vΦ, i])=(Piecewise[{{((((((Subscript[x, i-1])-x))^2) (-(Subscript[x, i-1])+3 (Subscript[x, i])-2 x))/((((Subscript[x, i])-(Subscript[x, i-1])))^3)), (Subscript[x, i-1])<x<(Subscript[x, i])}, {(((3 (Subscript[x, i])-(Subscript[x, i+1])-2 x) ((((Subscript[x, i+1])-x))^2))/((((Subscript[x, i])-(Subscript[x, i+1])))^3)), (Subscript[x, i])<x<(Subscript[x, i+1])}}])
 The slope shape function on the intervals [(Subscript[x, i-1]),(Subscript[x, i])] and [(Subscript[x, i]),(Subscript[x, i+1])] is 
 	(Subscript[sϕ, i])=(Piecewise[{{((((((Subscript[x, i-1])-x))^2) (x-(Subscript[x, i])))/((((Subscript[x, i])-(Subscript[x, i-1])))^2)), (Subscript[x, i-1])<x<(Subscript[x, i])}, {(((x-(Subscript[x, i])) ((((Subscript[x, i+1])-x))^2))/((((Subscript[x, i])-(Subscript[x, i+1])))^2)), (Subscript[x, i])<x<(Subscript[x, i+1])}}])

The second derivative of the value shape function on the intervals [(Subscript[x, i-1]),(Subscript[x, i])] and [(Subscript[x, i]),(Subscript[x, i+1])] is 
 	(Subscript[vΦ, i])''=(Piecewise[{{((6 (-(Subscript[x, i-1])-(Subscript[x, i])+2 x))/((((Subscript[x, i-1])-(Subscript[x, i])))^3)), (Subscript[x, i-1])<x<(Subscript[x, i])}, {((6 ((Subscript[x, i])+(Subscript[x, i+1])-2 x))/((((Subscript[x, i])-(Subscript[x, i+1])))^3)), (Subscript[x, i])<x<(Subscript[x, i+1])}}])
 The second derivative of the slope shape function on the intervals [(Subscript[x, i-1]),(Subscript[x, i])] and [(Subscript[x, i]),(Subscript[x, i+1])] is 
 	(Subscript[sΦ, i])''=(Piecewise[{{((-4 (Subscript[x, i-1])-2 (Subscript[x, i])+6 x)/((((Subscript[x, i-1])-(Subscript[x, i])))^2)), (Subscript[x, i-1])<x<(Subscript[x, i])}, {((-2 (Subscript[x, i])-4 (Subscript[x, i+1])+6 x)/((((Subscript[x, i])-(Subscript[x, i+1])))^2)), (Subscript[x, i])<x<(Subscript[x, i+1])}}])

We need to compute the following integrals for the RHS
	(∫_0^1)(Subscript[vϕ, i])[x]f[x]dx  and (∫_0^1)(Subscript[sϕ, i])[x]f[x]dx
We need to compute the following integrals for the LHS
	(Subscript[Avv, i,j])=(∫_0^1)(Subscript[vϕ, i])''[x](Subscript[vϕ, j])''[x]  and (Subscript[Ass, i,j])=(∫_0^1)(Subscript[sϕ, i])''[x](Subscript[sϕ, j])''[x]dx 
	and 
	(Subscript[Avs, i,j])=(∫_0^1)(Subscript[vϕ, i])''[x](Subscript[sϕ, j])''[x]dx
	and
	(Subscript[Asv, i,j])=(∫_0^1)(Subscript[sϕ, i])''[x](Subscript[vϕ, j])''[x]dx

we can precompute formulae for the LHS integrals if we are carefull.   As before, Avv and Ass are symmetric and the integrals are zero if |i-j|>1 because the shapefunctions do not overlap. We need to compute
	(Subscript[Avv, i,i]) and (Subscript[Avv, i,i+1])=(Subscript[Avv, i+1,i])
as well as
	 (Subscript[Ass, i,i]) and (Subscript[Ass, i,i+1])=(Subscript[Ass, i+1,i])

More Computation

First the building blocks for the diagonal entries.  Remember: the diagonal entries (excet the boundary ones) will have two pieces added together-one from each interval while the off-diagonal entries have a single piece from the overlap.

Clear[xp,xc,xn]
{Integrate[ 
vϕ[{xc,xp}]''[x]vϕ[{xc,xp}]''[x] ,
{x,xp,xc}],
Integrate[ 
sϕ[{xc,xp}]''[x]sϕ[{xc,xp}]''[x] ,
{x,xp,xc}],
Integrate[ 
sϕ[{xc,xp}]''[x]vϕ[{xc,xp}]''[x] ,
{x,xp,xc}]}

Next the building blocks for the off-diagonal entries

Clear[xp,xc,xn]
{Integrate[ 
vϕ[{xc,xp}]''[x]vϕ[{xp,xc}]''[x] ,
{x,xp,xc}],
Integrate[ 
sϕ[{xc,xp}]''[x]sϕ[{xp,xc}]''[x] ,
{x,xp,xc}],
Integrate[ 
sϕ[{xc,xp}]''[x]vϕ[{xp,xc}]''[x] ,
{x,xp,xc}]}

Book-keeping

We need to decide how we are going to arrange our pieces. I see two reasonable choices:

1) List of all the values and then all the slopes i.e
	(OverVector[u])={(Subscript[v, 1]),(Subscript[v, 2]),…,(Subscript[v, n]),(Subscript[s, 1]),(Subscript[s, 2]),…,(Subscript[s, n])}
2) list the points and slopes as they come
	(OverVector[u])={(Subscript[v, 1]),(Subscript[s, 1]),(Subscript[v, 2]),(Subscript[s, 2]),…,(Subscript[v, n]),(Subscript[s, n])}
as always sometimes one is easier than the other.

The first makes it easy to build the stiffness matrix as a block matrix.
	A=(({{Avv, Asv}, {Avs, Ass}}))
from the constituent pieces.