Numerical Analysis 7

Spring 2005

Problems with Bracketing Methods for Finding Roots of f(x)=0

The major problems are (1) if the function doesn't change sign, bracketing methods won't work, and (2) problems arise at multiple roots.

Plot[(x - 2)^2 (x + 2)^2, {x, -5, 5}];

Plot[x^5 - 1, {x, -3, 3}];

Error Considerations (Section 2.3)

Delta x and Delta y

Newton-Raphson and Secant Methods (Section 2.4)

Newton-Raphson

Outline

Assume that an initial approximation, (Subscript[p, 0]), is close to a root p of f(x)=0.  The next approximation, (Subscript[p, 1]), is the intersection of the tangent line to f(x) at the point ((Subscript[p, 0]), f((Subscript[p, 0]))) with the x-axis.  The equation of the tangent line is y-f((Subscript[p, 0]))=f '((Subscript[p, 0]))(x-(Subscript[p, 0])), and setting y=0 and solving for x gives x=(Subscript[p, 0])-((f((Subscript[p, 0])))/(f ' ((Subscript[p, 0])))).  So let (Subscript[p, 1])= (Subscript[p, 0])-((f((Subscript[p, 0])))/(f ' ((Subscript[p, 0])))), and in general (Subscript[p, k])= (Subscript[p, k - 1])-((f((Subscript[p, k - 1])))/(f ' ((Subscript[p, k - 1]))))

Theorem 2.5:  Assume that f ϵ (C^2)[a,b] and there is a number p ∈ [a,b] where f(p)=0.  If 
f '(p)≠0, then there is a δ>0 such that the sequence (Subscript[p, k]) defined by 
(Subscript[p, k]) = g((Subscript[p, k - 1]))=(Subscript[p, k - 1])-\(TraditionalForm\`\(f(p\_\(k - 1\))\)\/\(f\ '\ \((p\_\(k - 1\))\)\)\) for k=1,2,... will converge to p for any initial approximation (Subscript[p, 0])∈ [p-δ,p+δ].

An outline of the proof:  A Taylor Series expansion of f(x) about (Subscript[p, 0]) gives
f(x)=f((Subscript[p, 0])) + f ' ((Subscript[p, 0]))(x-(Subscript[p, 0]))+ ((f '' (c) (((x - (Subscript[p, 0])))^2))/(2!)) where c is between (Subscript[p, 0]) and x.  Substituting x=p into this equation and using the fact that f(p)=0 gives
0=f((Subscript[p, 0])) + f ' ((Subscript[p, 0]))(p-(Subscript[p, 0]))+ ((f '' (c) (((p - (Subscript[p, 0])))^2))/(2!))  
If (Subscript[p, 0]) is close enough to p the last term is much smaller than the sum of the first two terms.  Neglecting it gives
0≈f((Subscript[p, 0])) + f ' ((Subscript[p, 0]))(p-(Subscript[p, 0])), and solving this for p gives the iteration sequence.
Now  g(x)=x -((f(x))/(f ' (x))) so g ' (x) = 1 - \(TraditionalForm\`\(f\ '\ \((x)\) f\ '\ \((x)\) - \(f(x)\) f\ '' \((x)\)\)\/\((f\ '\ \((x)\))\)\^2\)= \(TraditionalForm\`\(\(f(x)\) f\ '' \((x)\)\)\/\((f\ '\ \((x)\))\)\^2\)  By hypothesis, f(p)=0 so g ' (p) = 0 and g ' (x) is continuous implies that there is a δδ>0 so that |g ' (x)|<1.  This means the fixed point iteration theorem is satisfied.

NewtonRaphson Computer Program

NewtonRaphson[f0_,x0_,delta_,epsilon_,max_,digits_,{leftx_,rightx_},{lowery_,uppery_}]
f0_ the function f(x),
x0_ the starting value,
delta_ estimate for the accuracy of the solution,
epsilon_ estimate for the magnitude of  f(x),
max_ the maximum number of iterations,
digits_ the number of digits displayed in the output
leftx_ is the leftmost x-value displayed on the graph
rightx_  is the rightmost x-value displayed on the graph
lowery_ is the minimum y-value displayed on the graph
uppery_ is the maximum y-value displayed on the graph
graphicstoggle_ is a toggle for graphs--they are shown if 1 and not shown if 0--unless you specify a zero, graphs are shown

NewtonRaphson[f0_, x0_, delta_, epsilon_, max_, 
        digits_, {leftx_, rightx_}, {lowery_, uppery_}, 
        graphicstoggle_:  1]:=
	Module[{f, k, small, cond, AbscissaInTolerance, 
            OrdinateInTolerance, RootInTolerance, Iterating},
		Set @@ {f[x_], f0};
		p0 = x0;
		Y0 = N[f[p0], digits];
		p1 = p0 + 1;
		k = 0;
	newp0 = p0;
	tablelist = {};
		cond = 0;
		small = N[(1/(10^20)), digits];
	gr1  = 
              Plot[f[x], {x, leftx, rightx}, 
                PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
                DisplayFunction -> Identity];
		While[cond == 0 && k < max,
			Df=N[f'[p0],digits];
			\(If[Df == 0, 
                    cond = DivByZero; 
				Dp = p1 - p0; 
				p1 = 
                      p0, 
				Dp = Y0\/Df; 
				p1 = p0 - Dp]\);
			Y1 = N[f[p1], digits];
		newp1 = p1;
gr2 = 
                    Show[Graphics[
                        Map[{PointSize[0.02], RGBColor[1, 0, 0], 
                              Point[#]} &, {{newp0, f[newp0]}, {newp1, 0}}]], 
                      DisplayFunction -> Identity];
gr3 = 
                    Show[Graphics[Line[{{newp0, f[newp0]}, {newp1, 0}}]], 
                      DisplayFunction -> Identity];
If[graphicstoggle == 1, 
                    Show[{gr1, gr2, gr3}, 
                      DisplayFunction -> $DisplayFunction]];
	tablelist = 
                    Append[tablelist, {k, newp0, newp1, Y1, 
                        Abs[newp0 - newp1]}];
			RelErr = 
                    N[Abs[Dp]/(Abs[p1] + small), digits];If[RelErr < delta, cond = AbscissaInTolerance];If[Abs[Y1] < epsilon, cond = OrdinateInTolerance];If[RelErr < delta && Abs[Y1] < epsilon, 
                    cond = RootInTolerance];p0 = p1;
			Y0 = Y1;
		newp0 = p1;
			k = k + 1;];
		If[cond == 0, cond = Iterating];Print["Sequence of Approximations for ", f[x]];
\(Print[
              TableForm[tablelist, 
                TableHeadings -> {None, {"k", "old", "new", "yval", "deltax"}}]]\);
tablelist
]

NewtonRaphson[x Sin[x] - 1, 1.7, 0.00001, 0.00001, 10, 
      6, {-2, 2}, {-5, 1}, 0]

Finding Square Roots Using Newton's Method

The function to use in the fixed point iteration procedure is g(x)=(1/2)(x+(N/x))

FixedPointIteration[(1/2) (3/x + x), 5, 0.01, 30, 
        12, {1, 5}, {0, 5}];

Or, a square root can be found by finding a zero of the function g(x)-x.

NewtonRaphson[(1/2) (x + (3/x)) - x, 10. , 0.00001, 0.00001, 10, 
      10, {0, 10}, {-5, 5}, 0]

Division by Zero Error

In the Newton-Raphson method formula, (Subscript[p, k])= (Subscript[p, k - 1])-\(TraditionalForm\`\(f(p\_\(k - 1\))\)\/\(f\ ' \((p\_\(k - 1\))\)\)\), the derivative at (Subscript[p, k - 1])occurs in the denominator.  If this derivative is zero the method cannot proceed.  One solution would be to take (Subscript[p, k - 1]) as the approximation to the root if f((Subscript[p, k - 1])) is close to zero.

Multiple Roots and Convergence

Definition: If f(x) and its derivatives \(TraditionalForm\`\(f\^'\)\)(x), \(TraditionalForm\`\(f\^''\)\)(x),...,(f^((M)))(x) are defined and continuous on an interval about x=p, we say that f(x)=0 has a root of order M at x=p if and only if 
f(p) = 0, \(TraditionalForm\`\(f\^'\)\)(p)=0, \(TraditionalForm\`\(f\^''\)\)(p)=0,...,(f^((M - 1)))(p)=0 but (f^((M)))(p)≠0 
A root of order 1 is called a simple root, one of order 2 is a double root, etc.

Example:

g[x_] := -8 + 20 x - 10 (x^2) - 13 (x^3) + 17 (x^4) - 7 (x^5) + 
        (x^6)

Plot[g[x], {x, -5, 5}];

D[g[x], {x, 3}] /. x -> 2

Lemma: If f(x)=0 has a root of order M at x=p, then there is a continuous function h(x) so that f(x)=(((x - p))^M)h(x) where h(p)≠0

Definition: Assume that (Subscript[p, n])-> p as n->∞ and set (Subscript[E, n])= p-(Subscript[p, n]) for n≥0.  If two positive constants, A≠0 and R>0 can be found such that 
(Underscript[(lim ), n -> ∞])\(TraditionalForm\`\(\(|\)\(p - p\_\(n + 1\)\)\(|\)\)\/\(\(|\)\(p - p\_n\)\( | \^R\)\)\)= (Underscript[(lim ), n -> ∞])\(TraditionalForm\`\(\(|\)\(E\_\(n + 1\)\)\(|\)\)\/\(\(|\)\(E\_n\)\( | \^R\)\)\)= A, then the sequence is said to converge to p with order of convergence R.  The number A is called the asymptotic error constant.

Consider the following examples, first considering the simple root at x=-2 and then the double root at x=1.

Example: f(x)=(x^3)-3x+2

Plot[(x^3) - 3  x + 2, {x, -5, 5}];

Simple Root at -2 with starting point -4.7

lll = 
      NewtonRaphson[(x^3) - 3  x + 2, -4.7, 0.01, 0.01, 20, 
        6, {-4, 0}, {-5, 5}, graphicstoggle_:  1]

appxs = Transpose[lll][[3]]

errs = Map[Abs[# - (-2)] &, appxs]

Map[(errs[[#]]/errs[[# - 1]]^2) &, 
      Drop[Range[Length[errs]], 1]]

Double Root at 1 with starting point 5.:

lll = 
      NewtonRaphson[(x^3) - 3  x + 2, 5. , 0.01, 0.01, 20, 
        6, {-4, 0}, {-5, 5}, graphicstoggle_:  1]

D[x^3 - 3  x + 2, x] /. x -> -2

D[x^3 - 3  x + 2, {x, 2}] /. x -> -2

appxs = Transpose[lll][[3]]

errs = Map[Abs[# - 1] &, appxs]

Map[(errs[[#]]/errs[[# - 1]]) &, 
      Drop[Range[Length[errs]], 1]]

Theorem 2.6: If the Newton-Raphson iteration produces a sequence (Subscript[p, n]) converging to root p of the function f(x) and p is a simple root, convergence is quadratic and 
|(Subscript[E, n + 1])|≈\(TraditionalForm\`\(\(|\)\(\(f\^''\)(p)\)\(|\)\)\/\(2 | \(f\^'\)(p) | \)\)|(Subscript[E, n])\( | \^2\) for n sufficiently large
If p is a multiple root of order M, convergence is linear and 
|(Subscript[E, n + 1])|≈((M - 1)/M)|(Subscript[E, n])| for n sufficiently large.

For the simple root at x=-2 of f(x)=(x^3) - 3  x + 2 shown above, ((|(Subscript[E, n + 1])|)/(|(Subscript[E, n])\( | \^2\)))≈\(TraditionalForm\`\(\(|\)\(\(f\^''\)(p)\)\(|\)\)\/\(2 | \(f\^'\)(p) | \)\)=(12/(2 (9)))= 2/3 and you can see this happening.  For the simple root at x=-1, also shown above, ((|(Subscript[E, n + 1])|)/(|(Subscript[E, n])|))≈((M - 1)/M)=((2 - 1)/2)= 1/2, and again, you see that this is the case from the computations done above.

Secant Method

The secant method finds the next approximation to a root by starting with 2 points close to the root.  The next point is found at the intersection of the line through these 2 points and the x-axis.  This process is repeated.

Clear[SecantMethod]

lll = 
      SecantMethod[(x^3) - 3  x + 2, -4.2, -4. , 0.001, 0.001, 20, 
        6, {-5, 0}, {-100, 3}]

\(\(SecantMethod[f0_, x0_, x1_, delta_, epsilon_, max_, 
          digits_, {leftx_, rightx_}, {lowery_, uppery_}, 
          graphicstoggle_:  1] := 
	
        Module[{f, k, small, cond, DivByZero, AbscissaInTolerance, 
            OrdinateInTolerance, RootInTolerance, Iterating, newp0, 
            newp1}, 
		Set @@ {f[x_], f0}; 
		p0 = x0; 
		Y0 = 
            N[f[p0], digits]; 
		p1 = x1; 
		Y1 = 
            N[f[p1], digits]; 
		k = 1; 
	
          newp0 = p0; 
	newp1 = p1; 
	
          tablelist = {{newp0, Y0}, {newp1, Y1}}; 
gr1\  = 
            Plot[f[x], {x, leftx, rightx}, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              DisplayFunction -> Identity]; 
		cond = 
            0; 
		small = N[1\/10\^20, digits]; 
	
          While[cond == 0 && 
              k < max, 

			Df = \(Y1 - Y0\)\/\(p1 - p0\); 
			If[Df == 0, cond = DivByZero; Dp = p1 - p0; p2 = p1, Dp = Y1\/Df; 
              p2 = p1 - Dp]; 
			Y2 = 
              N[f[p2], digits]; 
tablelist = 
              Append[tablelist, {p2, Y2}]; 
            gr2 = Show[
                Graphics[
                  Map[{PointSize[0.02], RGBColor[1, 0, 0], Point[#]} &, 
                    tablelist]], 
                DisplayFunction -> Identity]; 
gr3 = 
              Show[Graphics[Map[Line[#] &, addin[tablelist]]], 
                DisplayFunction -> Identity]; 
If[
              graphicstoggle == 1, 
              Show[{gr1, gr2, gr3}, 
                DisplayFunction -> $DisplayFunction]]; 

RelErr = 
              N[Abs[Dp]\/\(Abs[p2] + small\), digits]; 
            If[RelErr < delta, cond = AbscissaInTolerance]; 
            If[Abs[Y2] < epsilon, cond = OrdinateInTolerance]; 
            If[RelErr < delta && Abs[Y2] < epsilon, cond = RootInTolerance]; 
            p0 = p1; newp0 = p1; 
			Y0 = Y1; 
			p1 = p2; 
            newp1 = p2; 
			Y1 = Y2; 
			k = k + 1;]; 
		If[
            cond == 0, cond = Iterating]; 
          Print[" "]; 
		Print["Sequence of Approximations for ", 
            f[x]]; 
tablelist = 
            Transpose[
              Prepend[Transpose[tablelist], 
                  Table[m - 2, {m, 1, Length[tablelist]}]] /. \(-1\) -> 
                  0]; 
Print[
            TableForm[tablelist, 
              TableHeadings -> {None, {"Iteration", "pk", "f(pk)"}}]]
];\)\)

Functions used in this notebook

\(\(Bisection[f0_, a0_, b0_, delta_, 
          digits_, {leftx_, rightx_}, {lowery_, uppery_}, graphicstoggle_] := 
        Module[{a, b, c, Ya, Yb, Yc, f, k, cond, MaxIter, MethodDoesNotApply, 
            PerfectZero, RootInTolerance, tablelist, newa, newb, deltax, gr1, 
            gr2, gr3}, 
		Set @@ {f[x_], f0}; 
		a = 
            N[a0, digits]; 
		b = N[b0, digits]; 
	
          newa = a; newb = b; 
	
          tablelist = {}; 
		k = 0; 
		Ya = f[a]; 
		Yb = 
            f[b]; 
		dx = b - a; 
		cond = 0; 
If[
            Ya*Yb > 0, cond = MethodDoesNotApply; 
            Print["Starting Points Don't Have Opposite Signs"]; 
            Break[];]; 
gr1\  = 
            Plot[f[x], {x, leftx, rightx}, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              DisplayFunction -> Identity]; 
		MaxIter = 
            N[1 + Floor[\(Log[dx] - Log[delta]\)\/Log[2.0]]]; 
		While[
            cond == 0 && 
              k < MaxIter, 
		 (*If\ iteration\ count\ and\ dx\ are\ less\ than\ maximum\ continue*) 
		c = \(a + b\)\/2; 
			Yc = 
              f[c]; 

	 (*Add\ new\ point\ to\ list\ of\ points\ and\ display\ graphs*) 
	
            tablelist = 
              Append[tablelist, {k, newa, c, newb, Yc, 
                  deltax = Min[{Abs[newb - c], Abs[newa - c]}]}]; 
            gr2\  = Show[
                Graphics[
                  Map[{PointSize[0.03], RGBColor[1, 0, 0], Point[#]} &, {{c, 
                        0}}]], 
DisplayFunction -> 
                  Identity]; 
            gr3 = Show[
                Graphics[
                  Map[{PointSize[0.02], Point[#]} &, {{newa, 0}, {newb, 
                        0}}]], 
DisplayFunction -> 
                  Identity]; 
If[graphicstoggle == 1, 
              Show[{gr1, gr2, gr3}, 
                DisplayFunction -> $DisplayFunction]]; 
	
 (*\ 
              Check\ to\ see\ whether\ zero\ has\ been\ \(reached--\) 
                if\ so\ indicate\ a\ stop*) 
	
            If[Yc == 0, a = c; b = c; 
              cond = PerfectZero]; 

 (*\ 
              
              Check\ to\ see\ which\ point\ should\ be\ replaced\ by\ the\ new\ point\ and\ make\ the\ replacement*) 	
	
            If[Yb\ Yc > 0, b = c; newb = b; Yb = Yc, a = c; newa = c; 
              Ya = Yc]; 

 (*\ 
              Find\ the\ distance\ between\ endpoints\ a\ and\ b, \ 
              and\ increase\ iteration\ count\ by\ 1*) 
		\(dx = b - a; k = k + 1;\)]; 
	
          If[dx < delta && cond ≠ PerfectZero, 
            cond = RootInTolerance]; 
		Print["Sequence of Approximations for ", f[x]]; 
			Print[
            TableForm[tablelist, 
              TableHeadings -> {None, {"k", "left", "new midpoint", "right", "yval", "deltax"}}]];
];\)\)

\(FalsePosition[f0_, a0_, b0_, delta_, epsilon_, 
        digits_, {leftx_, rightx_}, {lowery_, uppery_}, graphicstoggle_:  1] := 
      Module[{a, b, c, Ya, Yb, Yc, f, k, cond, Iterating, RootInTolerance, 
          PerfectZero, MethodDoesNotApply, ptlist, nextpt, gr1, gr2, gr3, 
          gr4, gr5, gr6, gr7, newa, newb, 
          tablelist}, 
		Set @@ {f[x_], f0}; 
		a = 
          N[a0, digits]; 
		b = N[b0, digits]; 
		k = 0; newa = a; 
        newb = b; 
	
        tablelist = {}; 
		Ya = f[a]; 
		Yb = f[b]; 
		dx = 
          b - a; 
		cond = 0; 
		MaxIter = 
          100; 
		If[Ya\ Yb > 0, 
          cond = MethodDoesNotApply; 
          Print["Starting Points Don't Have Opposite Signs"]; 
          Break[];]; 
		gr1\  = 
          Plot[f[x], {x, leftx, rightx}, 
            PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
            DisplayFunction -> Identity]; 
		While[
          cond == 0 && 
            k < MaxIter, 
			dx = \(Yb\ \((b - a)\)\)\/\(Yb - Ya\); 
			c = b - dx; 
			Yc = f[c]; 
		tablelist = 
            Append[tablelist, {k, newa, c, newb, Yc, 
                Min[{Abs[newa - c], 
                    Abs[newb - 
                        c]}]}]; 
			 (*\(Print["c"\_k, " = ", 
                N[c, digits], ",   |dc| ≤ ", 
                N[dx, digits]];\)*) 
			If[Yc == 0, a = c; b = c; 
            cond = PerfectZero]; 
gr2 = 
            Show[Graphics[
                Map[{PointSize[0.02], Point[#]} &, {{newa, Ya}, {newb, 
                      Yb}}]], 
              DisplayFunction -> Identity]; 
gr3 = 
            Show[Graphics[Line[{{newa, Ya}, {newb, Yb}}]], 
              DisplayFunction -> Identity]; 
If[
            Yb\ Yc > 0, b = c; newb = b; Yb = Yc, a = c; newa = a; 
            Ya = Yc]; 
	
          gr4 = Show[
              Graphics[
                Map[{PointSize[0.03], RGBColor[1, 0, 0], Point[#]} &, {{c, 
                      0}}]], DisplayFunction -> Identity]; 
          If[graphicstoggle == 1, 
            Show[{gr1, gr2, gr3, gr4}, 
              DisplayFunction -> $DisplayFunction]]; 
If[Abs[dx] < delta && Abs[Yc] < epsilon, 
            cond = RootInTolerance]; k = k + 1;]; 
        If[cond == 0, 
          cond = Iterating]; 
		Print["Sequence of Approximations for ", f[x]]; 
Print[
          TableForm[tablelist, 
            TableHeadings -> {None, {"k", "left", "newpoint", "right", "yval", "deltax"}}]]
]\)

\(\(FixedPointIteration[g0_, X0_, tolerance_, max_, 
          digits_, {leftx_, rightx_}, {lowery_, uppery_}] := 
	
        Module[{g, k, x0, x1, relerr, plist, doubleplist, gr1, gr2, gr3, 
            gr4}, 
 (*Initial\ Values\ for\ Variables*) 
	
          Set @@ {g[x_], g0}; 
		x0 = 
            N[X0, digits]; 
	\ \ plist = {x0}; 
		k = 
            0; 
		err = 1. ; 
		relerr = 1; 
		small = 
            N[1\/10\^20, 
              digits]; 
Print[\*""Iterates (Subscript[p, k]),  "", "    Relative Error |Δp|"]; 
		Print["p"\_0, " = ", N[x0, digits]]; 
 (*Iteration*) 
		While[
            k < max && 
              tolerance ≤ 
                relerr, 
			Module[{}, 
				x1 = 
                g[x0]; 
				k = k + 1; 
				err = 
                Abs[x1 - x0]; 
				Print["p"\_k, " = ", 
                N[x1, digits], ",   |Δp| ≤ ", 
                N[err, digits]]; 
			relerr = \(2  Abs[x1 - x0]\)\/\(Abs[x1] + small\); 
				x0 = x1; 
              plist = Append[plist, x1]]]; 
		p = 
            x0; 
		Print[" "]; 
 (*Collect\ approximations\ into\ a\ list\ to\ form\ points\ for\ graphing*) 
	
          doubleplist = 
            Flatten[Transpose[{plist, plist}]]; 
	
          pts = Transpose[{Drop[doubleplist, {Length[doubleplist]}], 
                Prepend[Drop[doubleplist, 2], 
                  0]}]; 
 (*Print\ Results*) 
		Print["The function is g[x] = ", g[x]]; 
          Print["Approximation to a fixed point of g[x] is:"]; 
          Print["  p  = ", N[p, digits]]; 
          Print["g[p] at the approximated fixed point= ", 
            N[g[p], digits]]; 
          Print["|Δp| ≤ ", 
            N[err, digits]]; 
 (*Graph\ Results*) 
	
          gr1\  = Plot[g[x], {x, leftx, rightx}, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              PlotStyle -> RGBColor[1, 0, 0], 
              DisplayFunction -> Identity]; 
		gr2\  = 
            Plot[x, {x, leftx, rightx}, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              PlotStyle -> RGBColor[0, 1, 0], 
              DisplayFunction -> Identity]; 
	
          gr3\  = ListPlot[pts, 
              PlotRange -> {{leftx, rightx}, {lowery, uppery}}, 
              PlotStyle -> {PointSize[0.01]}, 
              DisplayFunction -> Identity]; 
	
          gr4 = Show[Graphics[Line[pts]], 
              DisplayFunction -> Identity]; 
		Show[{gr1, gr2, gr3, 
              gr4}, DisplayFunction -> $DisplayFunction];];\)\)

addin[
        lll_] := 
  
  
Module[{temp, rest, newones}, 
temp = 
          Take[lll, 2]; 
rest = 
          Drop[lll, 2]; 
newones = 
          Transpose[rest]; 
newones = 
          Append[newones, 
            Table[0, {Length[
                  newones[[1]]]}]]; 
Partition[
          Partition[
            Flatten[Prepend[
                Transpose[{newones[[1]], newones[[3]], 
                    newones[[1]], newones[[2]]}], temp]], 2], 2, 
          1]
]