Off[General::spell] ;       Off[
      General::spell1]  ;      Off[General::luc];

(*        Rechengenauigkeit   vorgeben  ipl = 1, 2, 3 ,4  oder 5    *)

ipl =  1;  
(*     Nummern    
            der  auszugebenden Kurven festlegen       *)
If [ 
        ipl  ==  1,  { ipl1  =  1;  
          ipl2  =  1;   ipl3  =  1;   ipl4  =  1, 
          ipl5  =  1 }];
If [ 
        ipl  ==  2,  { ipl1  =  1;  
          ipl2  =  2;   ipl3  =  1;   ipl4  =  1, 
          ipl5  =  1 }];
If [ 
        ipl  ==  3,  { ipl1  =  1;  
          ipl2  =  2;   ipl3  =  3;   ipl4  =  1, 
          ipl5  =  1 }];
If [ 
        ipl  ==  4,  { ipl1  =  1;  
          ipl2  =  2;   ipl3  =  3;   ipl4  =  4, 
          ipl5  =  1 }];
If [ 
        ipl  ==  5,  { ipl1  =  1;  
          ipl2  =  2;   ipl3  =  3;   ipl4  =  4, 
          ipl5  =  5 }];

listgen= { 8,12,16,24 ,32};        
   $MaxPrecision = gen =  listgen[[ipl]];
   $MinPrecision = gen;

(*            Approximation ( Methode der kleinsten Quadrate nach Gauß                       *)
(*                       Approximationsfunktionen  ( Polynome vom Grade  m )                              *)
(*                       für n+1 vorgegene Knotenpunkte  mit Hilfe der  Normalgleichungen              *)

n = 20;

SeedRandom[1243]
f5[x_] := 0.01  (x^5)-0.525  (x^4)+10.41 
           (x^3)-96.47  (x^2)+411.17  
          x-530.26+10*
          Random[]-5

Do[{xp[j] = SetPrecision[j, gen],  yp[j] = f5[j]}, {j, 0, 
        n}]

<< Graphics`Colors`
liste1 = {Red, HotPink, Green, Apricot, Brown, DarkGreen, Cobalt, 
          Brick, Orange, Magenta, IndianRed, ForestGreen, Red, HotPink, 
          Green, Apricot, Brown, DarkGreen, Cobalt, Brick, 
          Orange};

\(linienplot = 
      ListPlot[Table[{xp[j], yp[j]}, {j, 0, n}], PlotJoined\  -> \ True, 	
        PlotRange -> {{0, 21}, {\(-210\), 210}}, 
	PlotStyle -> Blue, 
        AspectRatio -> 0.6, AxesLabel -> {"> x", "  ^ y"}]\)

\(punktplot = 
      ListPlot[Table[{xp[j], yp[j]}, {j, 0, n}], PlotJoined\  -> \ False, 	
        PlotRange -> {{0, 21}, {\(-210\), 210}}, 
	PlotStyle -> Blue, 
        Prolog\  -> \ AbsolutePointSize[5], AspectRatio -> 0.6, 
        AxesLabel -> {"> x", "  ^ y"}]\)

Show[linienplot, punktplot, 
      Prolog  ->  AbsolutePointSize[5]]

m  =  20;

\(\(\(Do[{A[j, 0] = SetPrecision[1, gen], 
        Do[A[j, i] = SetPrecision[xp[j]^i, gen], {i, 1, m}]}, {j, 0, n}]

    \(Do[Amat[k]\  = \ Table[A[j, i], {j, 0, n}, {i, 0, k}], {k, 1, 
          m}];\)\)\(\ \)\)\)

(*   Amat[5]  //  MatrixForm  ,    
          Kontrollausgabe   *)

Do[ATmat[k] = Transpose[Amat[k]], {k, 1, m}];

(*  ATmat[5]  //  MatrixForm ,   Kontrollausgabe   *)

yvek =  
        Table[SetPrecision[yp[j], gen], {j, 0, n}];

(*  yvek  //  MatrixForm   ,    
          Kontrollausgabe   *)

Do[Nmat[k]  =  SetPrecision[ATmat[k] . Amat[k], gen], {k, 1, 
        m}]

(*  Nmat[5]  //  MatrixForm  ,    
            Kontrollausgabe   *)

Do[bvek[k]  = SetPrecision[ATmat[k] . yvek, gen], {k, 1, 
        m}]

(*  bvek[5]  //  MatrixForm    ,    
            Kontrollausgabe   *)

Do[
        Ninv[k] = SetPrecision[Inverse[Nmat[k]], gen], {k, 1, 
          m}];

(*  Ninv[5]  //  MatrixForm ,   Kontrollausgabe    *)

Do[cvek[k] = SetPrecision[Ninv[k] . bvek[k], gen], {k, 1, 
        m}]

(*  cvek[9]  //  MatrixForm   ,    Kontrollausgabe   *)

poly[kk_, 
        xx_] := (∑_(i = 1)^kk) cvek[kk][[i + 1]]*xx^i + cvek[
            kk][[1]]

Do[{Dquad[ipl, 
                k] = (∑_(j = 0)^n)(yp[j] - poly[k, xp[j]])^2,Print["Ordnung der Ersatzfunktion: ", 
              k, ", Abweichungsquadrate: ", 
              ScientificForm[Dquad[ipl, k]], "   Gauss"],Kurv[k]=Plot[ poly[k, xx], {xx, 0, 21}, 
                PlotRange -> {{0, 21}, {-210, 210}}, 
                AspectRatio -> 
                  0.6, 
PlotPoints -> 40, 
                AxesLabel -> { "-> X", "  ^ Y"},  
                PlotStyle -> liste1[[k]]]},{k, 1, m}]

Do[{Print["Ordnung der Ersatzfunktion: ", 
          k, ", Abweichungsquadrate: ", 
          ScientificForm[Dquad[ipl, k]], "   Gauss"], 
        Show[Kurv[k], linienplot, punktplot, 
          Prolog  ->  AbsolutePointSize[5]]}, {k, 1, 20}]

(*  für l  =  
            1,  . . . ,  m    *)

(*     Vorwärtsrechnung ( 
                erste i - Schleife   ) ,  Rückwärtsrechnung ( 
                zweite i - Schleife )    *)

Do[{nl = l + 1,
      	Do[
              Do[ Achol[i, j] = Nmat[l][[i, j]], {i, 1, nl}], {j, 1, 
                nl}],
		\(Do[{Do[{sumAjiAjk = SetPrecision[0. , gen], 
                    Do[sumAjiAjk = sumAjiAjk + Achol[j, i]*Achol[j, k], {j, 
                        1, i - 
                          1}], 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	Achol[
                        i, k] = \((Achol[i, k] - sumAjiAjk)\)/
                        Achol[i, i]}, {i, 1, k - 1}], 
                sumAjk2 = 
                  SetPrecision[0. , 
                    gen], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Do[
                  sumAjk2 = sumAjk2 + Achol[j, k]*Achol[j, k], {j, 1, 
                    k - 1}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AkkminsumAjk2 = SetPrecision[Achol[k, k] - sumAjk2, gen], 
                Achol[k, k] = 
                  SetPrecision[Sqrt[AkkminsumAjk2], 
                    gen]}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {k, 1, nl}]\),		
	Do[{sumAb = SetPrecision[0. , gen], 
                zwb[i] = SetPrecision[bvek[l][[i]], gen],  
                Do[sumAb = sumAb + Achol[j, i]*zwb[j], {j, 1, i - 1}], 
                zwb[i] = 
                  SetPrecision[(zwb[i] - sumAb)/Achol[i, i], gen]}, {i, 
                1, nl}],
	 \(Do[{sumAX = SetPrecision[0. , gen], 
                Do[sumAX = sumAX + Achol[i, j]*zwc[j], {j, i + 1, 
                    nl}], 
	\ \ zwc[i] = 
                  SetPrecision[\((zwb[i] - sumAX)\)/Achol[i, i], gen]}, {i, 
                nl, 1, \(-1\)}]\) ,  cvec[l] = Table[zwc[i], {i, 1, nl}],
       
    If[
              AkkminsumAjk2  ≤  0, { 
                Print[" Die Matrix ist bei der Ordnung n = ", 
                  nl, " nicht mehr positiv definit !!!  "]}]},    
         {l, 1, m}]

polychol[kk_, 
        xx_] := (∑_(i = 1)^kk) cvec[kk][[i + 1]]*xx^i + cvec[
            kk][[1]]

Do[{Dquadchol[ipl, 
                  k] = (∑_(j = 0)^n)(yp[j] - 
                        polychol[k, xp[j]])^2,Print["Ordnung der Ersatzfunktion: ", 
                k, ", Abweichungsquadrate: ", 
                ScientificForm[Dquadchol[ipl, k]], "   Cholesky"],Kurvchol[k]=Plot[ polychol[k, xx], {xx, 0, 21}, 
                  PlotRange -> {{0, 21}, {-210, 210}}, 
                  AspectRatio -> 
                    0.6, 
PlotPoints -> 40, 
                  AxesLabel -> { "-> X", "  ^ Y"},  
                  PlotStyle -> liste1[[k]]]},{k, 1, m}]

Do[{Print["Ordnung der Ersatzfunktion: ", 
          k, ", Abweichungsquadrate: ", 
          ScientificForm[Dquadchol[ipl, k]],  "   Cholesky"], 
        Show[Kurvchol[k], linienplot, punktplot, 
          Prolog  ->  AbsolutePointSize[5]]}, {k, 1, m}]

liste2 = {Red, HotPink, Green, DarkGreen, Cobalt, Blue, Brick, Brown, 
          Orange, Magenta, Apricot, IndianRed, ForestGreen, };

Print["                   Logarithmus der Abweichungsquadrate bei Gau\
ß"]
\(quadlinien[ipl] = 
      ListPlot[Table[{k, Log[10, Dquad[ipl, k]]}, {k, 1, m}], 
        PlotJoined\  -> \ True, 	PlotRange -> {{0, 21}, {\(-1\), 11}}, 
	
        PlotStyle -> liste2[\([2  ipl - 1]\)], AspectRatio -> 0.7, 
        PlotLabel -> " ", 
        AxesLabel -> {"> k", "  ^ y"}]\)

Print["                   Logarithmus der Abweichungsquadrate bei \
Gauß"];
\(\(quadpunkte[ipl] = 
        ListPlot[Table[{k, Log[10, Dquad[ipl, k]]}, {k, 1, m}], 
          PlotJoined\  -> \ False, 	
          PlotRange -> {{0, 21}, {\(-1\), 11}}, 
	
          PlotStyle -> liste2[\([2  ipl - 1]\)], 
          Prolog\  -> \ AbsolutePointSize[5], AspectRatio -> 0.7, 
          PlotLabel -> " ", 
          AxesLabel -> {"> k", "  ^ y"}];\)\)

Print["                   Logarithmus der Abweichungsquadrate bei \
Cholesky"];
\(\(quadlinienchol[ipl] = 
        ListPlot[Table[{k, Log[10, Dquadchol[ipl, k]]}, {k, 1, m}], 
          PlotJoined\  -> \ True, 	PlotRange -> {{0, 21}, {\(-1\), 11}}, 
	
          PlotStyle -> liste2[\([2*ipl]\)], AspectRatio -> 0.7, 
          PlotLabel -> " ", 
          AxesLabel -> {"> k", "  ^ y"}];\)\)

\(\(quadpunktechol[ipl] = 
        ListPlot[Table[{k, Log[10, Dquadchol[ipl, k]]}, {k, 1, m}], 
          PlotJoined\  -> \ False, 	
          PlotRange -> {{0, 21}, {\(-1\), 11}}, 
	
          PlotStyle -> liste2[\([2  ipl]\)], 
          Prolog\  -> \ AbsolutePointSize[5], AspectRatio -> 0.7, 
          PlotLabel -> "Abweichungsquadrate , Cholesky", 
          AxesLabel -> {"> k", "  ^ y"}];\)\)

Show[quadlinien[ipl1], quadpunkte[ipl1], quadlinienchol[ipl1], 
        quadpunktechol[ipl1], quadlinien[ipl2], quadpunkte[ipl2], 
        quadlinienchol[ipl2], quadpunktechol[ipl2], quadlinien[ipl3], 
        quadpunkte[ipl3], quadlinienchol[ipl3], quadpunktechol[ipl3], 
        quadlinien[ipl4], quadpunkte[ipl4], quadlinienchol[ipl4], 
        quadpunktechol[ipl4], 
quadlinien[ipl5], 
        quadpunkte[ipl5], quadlinienchol[ipl5], 
        quadpunktechol[ipl5], 
Prolog  ->  
          AbsolutePointSize[5]];

\(\(If[
        ipl == \ 
          5, {Print["                   Logarithmus der Abweichungsquadrate bei Gauß"]; 
quadlinien[ipl] = 
            ListPlot[Table[{k, Log[10, Dquad[ipl, k]]}, {k, 1, m}], 
              PlotJoined\  -> \ True, 	
              PlotRange -> {{0, 21}, {\(-21\), 7}}, 
	
              PlotStyle -> liste2[\([2*ipl - 1]\)], AspectRatio -> 0.7, 
              PlotLabel -> " ", 
              AxesLabel -> {"> k", "  ^ y"}]}];\)\)

If[
        ipl ==  
          5, {Print["                   Logarithmus der Abweichungsquadrate \
bei Cholesky"]; 
          quadpunkte[ipl] = 
            ListPlot[Table[{k, Log[10, Dquad[ipl, k]]}, {k, 1, m}], 
              PlotJoined  ->  False, 
              PlotRange -> {{0, 21}, {-21, 
                    7}}, 
PlotStyle -> 
                liste2[[2*ipl - 1]], Prolog  ->  AbsolutePointSize[5], 
              AspectRatio -> 0.7, PlotLabel -> " ", 
              AxesLabel -> {"> k", "  ^ y"}]}];

\(\(If[
        ipl == \ 
          5, {Print["                   Logarithmus der Abweichungsquadrate bei Cholesky"]; 
quadlinienchol[pil] = 
            ListPlot[Table[{k, Log[10, Dquadchol[ipl, k]]}, {k, 1, m}], 
              PlotJoined\  -> \ True, PlotRange -> {{0, 21}, {\(-21\), 7}}, 

              PlotStyle -> liste2[\([2*ipl]\)], AspectRatio -> 0.7, 
              PlotLabel -> " ", 
              AxesLabel -> {"> k", "  ^ y"}]}];\)\)

\(\(If[ipl == \ 5, 
        quadpunktechol[ipl] = 
          ListPlot[Table[{k, Log[10, Dquadchol[ipl, k]]}, {k, 1, m}], 
            PlotJoined\  -> \ False, 	
            PlotRange -> {{0, 21}, {\(-21\), 7}}, 
	
            PlotStyle -> liste2[\([2*ipl]\)], 
            Prolog\  -> \ AbsolutePointSize[5], AspectRatio -> 0.7, 
            PlotLabel -> "Abweichungsquadrate , Cholesky", 
            AxesLabel -> {"> k", "  ^ y"}]];\)\)

If[ipl ==  5, 
        Show[quadlinien[ipl5], quadpunkte[ipl5], quadlinienchol[ipl5], 
          quadpunktechol[ipl5], 
          Prolog  ->  AbsolutePointSize[5]]];