(*     Interpolation:                                                                                                         *)
(*                 Parametrische Darstellung   x= φ(x) und y = ψ(x) )                      *)                                         
(*                       ( für geschlossene Kurven )                                                  *)
(*                              1.  Klassische  Interpolation   ( Newton )                       *)
(*                              2.  Natürliche kubische Spline-Funktion                      *)

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

lauf = 1;                         (*  Lauf mit der Nummer lauf  ,  
          alle vorhergehenden Läufe bleiben erhalten ! *)

sch =  "V";

If[sch == "N", 
      If[index ≥ 0 And index ≤ 
          m  ,    
TableForm[
          Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
              PaddedForm[xp[i], {3, 2}], 
"    yp[", 
              ind, "] = ", PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
          TableSpacing -> {1, 0}]]]

(*  Stützpunkte 0,  .. , m  eingeben        *)

\(\(If[\ 
        sch == "V", {\ \ m = 
            20, 
xp[0] = \(-4. \); \ \ \ yp[0] = 
            0. ; \ \ \ xp[1] = \(-3.95\); \ \ \ \ \ yp[1] = 
            0.5; \ \ \ \ \ \ xp[2] = \(-3.5\); \ \ \ \ \ \ yp[2] = 
            1.5; \ \ \ \ \ xp[3] = \(-3. \); \ \ \ yp[3] = 
            2. ; \ \ \ xp[4] = \ \(-2\); \ \ \ \ \ \ \ \ \ \ yp[4] = 
            2.6; \ \ \ \ \ \ xp[5] = \ \(-1\); \ \ \ \ \ \ \ \ \ yp[5] = 
            2.9; \ \ \ \ \ \ \ \ xp[6] = \ \ \ 0. ; \ \ \ yp[6] = 
            3. ; \ \ xp[7] = \ \ \ 2; \ \ \ \ \ \ \ \ \ \ \ yp[7] = 
            2.6; \ \ \ \ \ \ \ xp[8] = \(+3. \); \ \ \ \ \ \ \ \ yp[8] = 
            2. ; 
xp[9] = \ 3.75; \ \ \ \ \ \ yp[9] = 
            1; \ \ \ \ \ \ \ \ \ \ xp[10] = \(+4. \); \ \ \ \ \ \ \ yp[10] = 
            0. ; \ \ \ \ \ \ \ \ xp[11] = 
            3.5; \ \ \ \ \ yp[11] = \(-1.5\); \ \ \ \ \ \ 

          xp[12] = \(+3. \); \ \ \ \ \ \ yp[12] = \(-2. \); \ \ \ \ xp[
              13] = \ \ \ 1; \ \ \ \ \ \ \ \ \ yp[13] = \(-2.9\); \ \ xp[
              14] = \ \ \ 0. ; \ \ \ yp[14] = \(-3. \); 
xp[
              15] = \ \(-2\); \ \ \ \ \ \ \ \ yp[15] = \(-2.6\); \ 
          xp[16] = \(-3. \); \ \ \ \ \ yp[16] = \(-2. \); \ \ \ \ \ xp[
              17] = \(-3.5\); \ \ \ yp[
              17] = \(-1.5\); \ \ \ \ \ \ 
xp[
              18] = \(-3.75\); \ 
          yp[18] = \(-1. \); \ \ xp[19] = \(-3.95\); \ \ \ \ yp[
              19] = \(-0.5\); \ \ xp[20] = \(-4. \); \ \ \ yp[20] = 
            0. ;}];\)\)

(*  Vorgegebene Stützpunkte 0,  .. , m       *)

If[sch == "V", 
      TableForm[
        Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
            PaddedForm[xp[i], {3, 2}], "    yp[", ind, "] = ", 
            PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
        TableSpacing -> {1, 0}]]

index  = 9 ;                  
  \(If[sch == "L", 
          If[index ≥ 0\ And\ index ≤ 
              m\ , {m = 
                m - 1; 
\ \ \ \ \ \ \ \ \ Do[{\ xp[j] = xp[j + 1], 
                  yp[j] = yp[j + 1]}, {j, index, m}]}]]\);

(*   Stützpunkte   0,..,m ,  nach dem Entfernen von Punkten    *)

If[sch == "L", 
      If[index ≥ 0 And index ≤ 
          m  ,    
TableForm[
          Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
              PaddedForm[xp[i], {3, 2}], 
"    yp[", 
              ind, "] = ", PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
          TableSpacing -> {1, 0}]]]

index  =   1 ;                  
If[sch == "E",If[index ≥ 0 And index ≤ 
                    m,{m = m + 1;Do[{xp[j], yp[j]}, {j, 0, index - 1}];
   Do[{xp[j] = xp[j - 1], yp[j] = yp[j - 1]}, {j, 
                          m, index + 1, -1}];
       xp[index] =   - 3.75                                         ;
        yp[index] =    1.0                                       }]];

nummer[lauf] = m;

(*   Stützpunkte   0,..,m ,  nach dem Einfügen von Punkten    *)

If[sch == "E", 
      If[index ≥ 0 And index ≤ 
          m,    
TableForm[
          Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
              PaddedForm[xp[i], {3, 2}], 
"    yp[", 
              ind, "] = ", PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
          TableSpacing -> {1, 0}]]]

index  =   19 ;                  
    If[sch == "K",If[index ≥ 0 And index ≤ 
                    m,{       xp[index] =    -3.90                                         ;
        yp[index] =    -0.5                                       }]];

(*   Stützpunkte   0,..,m ,  nach dem Einfügen von Punkten    *)

If[sch == "K", 
      If[index ≥ 0 And index ≤ 
          m,    
TableForm[
          Table[{ind = PaddedForm[i, 2], "     xp[", ind, "] = ", 
              PaddedForm[xp[i], {3, 2}], 
"    yp[", 
              ind, "] = ", PaddedForm[yp[i], {3, 2}]}, {i, 0, m}], 
          TableSpacing -> {1, 0}]]]

tp[0] = 0;

    Do[{Delt[k - 1] = 
          Sqrt[(xp[k] - xp[k - 1])^2 + (yp[k] - yp[k - 1])^2], 
        tp[k] = tp[k - 1] + Delt[k - 1]}, {k, 1, m}]

nd = 400;

Do[{tj[j] = tp[0] + j*(tp[m] - tp[0])/nd,  
        tel[j] = 
          8. *ArcTan[1]/nd*
            j, 
                                          xel[j] = - 4*Cos[tel[j]], 
        yel[j] = 3*Sin[tel[j]]}, {j, 0, nd}]

<< Graphics`Colors`

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

\(ellipsplot = 
      ListPlot[Table[{xel[j], yel[j]}, {j, 0, nd}], 
        PlotJoined\  -> \ True, 
	
        PlotRange -> {{\(-5\), 5}, {\(-4\), 4}}, 	PlotStyle -> Brown, 
        AspectRatio -> 0.5, AxesLabel -> {"> x", "  ^ y"}]\)

\(punktplot = 
      ListPlot[Table[{xp[k], yp[k]}, {k, 0, m}], PlotJoined\  -> \ False, 	
        PlotRange -> {{\(-5\), 5}, {\(-4\), 4}}, 
	PlotStyle -> Blue, 
        Prolog\  -> \ AbsolutePointSize[5], AspectRatio -> 0.5, 
        AxesLabel -> {"> x", "  ^ y"}]\)

\(linienplot = 
      ListPlot[Table[{xp[k], yp[k]}, {k, 0, m}], PlotJoined\  -> \ True, 	
        PlotRange -> {{\(-5\), 5}, {\(-4\), 4}}, 
	PlotStyle -> Blue, 
        Prolog\  -> \ AbsolutePointSize[5], AspectRatio -> 0.5, 
        AxesLabel -> {"> x", "  ^ y"}]\)

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

FunkDivDiff[xpform_, 
        ypform_] :=  (
Do[
          DivDiffret[k, 
              1] = (ypform[k + 1] - ypform[k])/(xpform[k + 1] - 
                  xpform[k]), {k, 0, m - 1}]; 
Do[
          Do[DivDiffret[k, 
                j] = (DivDiffret[k + 1, j - 1] - 
                    DivDiffret[k, j - 1])/(xpform[k + j] - 
                    xpform[k]), {k, 0, m - j}], {j, 2, m}]; 
        DivDiffmatret = 
          Table[DivDiffret[0, j], {j, 1, m}]; {DivDiffmatret})

DivDiffxmat = FunkDivDiff[tp, xp] ;

DivDiffx = DivDiffxmat[[1]];

DivDiffymat = FunkDivDiff[tp, yp] ;

DivDiffy = DivDiffymat[[1]];

\(Do[{xnew[j] = xp[0], 
        prox = tj[j] - 
            tp[0], 
Do[{xnew[j] = 
              xnew[j] + prox*DivDiffx[\([i]\)], 
            prox = prox*\((tj[j] - tp[i])\)}, {i, 1, 
            m}], 
ynew[j] = yp[0], 
        proy = tj[j] - tp[0], 
	
        Do[{ynew[j] = ynew[j] + proy*DivDiffy[\([i]\)], 
            proy = proy*\((tj[j] - tp[i])\)}, {i, 1, m}]}, {j, 0, 
        nd}]\)

\(newtonplot = 
      ListPlot[Table[{xnew[j], ynew[j]}, {j, 0, nd}], 
        PlotJoined\  -> \ True, 
	
        PlotRange -> {{\(-5\), 5}, {\(-4\), 4}}, 	PlotStyle -> Red, 
        AspectRatio -> 0.5, AxesLabel -> {"> x", "  ^ y"}]\)

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

(*           Natürliche Spline - 
            Interpolation                         *)

FunkNatSpl[xpform_, ypform_]:= (
                                        
                                               
Du[0] = xpform[1] - xpform[0];
    Do[{Du[k] = xpform[k + 1] - xpform[k], 
              Dh[k] = 2 (Du[k - 1] + Du[k])}, {k, 1, m - 1}];
                                        
Ch[1] = Sqrt[Dh[1]];
    Do[{Cn[k - 1] = Du[k - 1]/Ch[k - 1], 
              Ch[k] = Sqrt[Dh[k] - Cn[k - 1]^2]}, {k, 2, m - 1}];
                      
                      
Dv[0] = ypform[1] - ypform[0];       Dv[1] = ypform[2] - ypform[1];
    Dr[1] = 3 (Dv[1]/Du[1] - Dv[0]/Du[0]);
    Z[1] = Dr[1]/Ch[1];
    \(Do[{Dv[k] = ypform[k + 1] - ypform[k], 
		Dr[k] = 
                3 \((Dv[k]/Du[k] - Dv[k - 1]/Du[k - 1])\), 
		Z[
                  k] = \((Dr[k] - Z[k - 1]*Cn[k - 1])\)/Ch[k]}, 
\ \ \ {k, 
              2, m - 1}]\);
                           
                           
B[m] = 0;  B[m - 1] = Z[m - 1]/Ch[m - 1];
    Do[
            B[k] = (Z[k] - B[k + 1]*Cn[k])/Ch[k], {k, m - 2, 
              1, -1}]; B[0] = 0;
                          
\(Do[{A[
                  k] = \((B[k + 1] - B[k])\)/\((3 
                       Du[k])\), 
		\ \ \ \ \ \ \ \ \ \ \ \ \ Cc[k] = 
                Dv[k]/Du[k] - Du[k]*\((B[k + 1] + 2  B[k])\)/3, \ 
              Dc[k] = ypform[k]}, {k, 0, m - 1}]\);
  Aret = Table[A[k],  {k, 0, m - 1}];      Bret = Table[B[k], {k, 0, m}];   
  Cret = Table[Cc[k], {k, 0, m - 1}];    Dret = Table[Dc[k], {k, 0, m - 1}];
{Aret, Bret, Cret, Dret})

ABCDxmat = FunkNatSpl[tp, xp]  ;

Anatx = ABCDxmat[[1]] ;

Bnatx = ABCDxmat[[2]];                       (* Koeffizienen Bk  f  ür die x-Werte *)

Cnatx = ABCDxmat[[3]];                        (* Koeffizienen Ck  f  ür die x-Werte *)

Dnatx = ABCDxmat[[4]];                        (* Koeffizienen Dk  f  ür die x-Werte *)

ABCDymat = FunkNatSpl[tp, yp] ;

Anaty = ABCDymat[[1]];

Bnaty = ABCDymat[[2]];                       (* Koeffizienen Bk  f  ür die x-Werte *)

Cnaty = ABCDymat[[3]];                        (* Koeffizienen Ck  f  ür die x-Werte *)

Dnaty = ABCDymat[[4]] ;                       (* Koeffizienen Dk  f  ür die x-Werte *)

\(Do[{tint = tp[0], 
	
        Do[{Dnt[k] = tp[k + 1] - tp[k], 
            If[\((tj[j] ≥ 
                    tint\ )\)\  ∧ \ \ \((tj[j] ≤ 
                    tint + Dnt[k]\ )\), {knt = k, Break[]}\ , 
              tint = tint + Dnt[k]]}, {k, 0, m - 1}], 
	
        Dntmj = tj[j] - tp[knt], 
\ \ \ \ \ xnat[j] = 
          Anatx[\([knt + 1]\)]*Dntmj^3 + Bnatx[\([knt + 1]\)]*Dntmj^2 + 
            Cnatx[\([knt + 1]\)]*Dntmj + Dnatx[\([knt + 1]\)]\ , 
	
        ynat[j] = 
          Anaty[\([knt + 1]\)]*Dntmj^3 + Bnaty[\([knt + 1]\)]*Dntmj^2 + 
            Cnaty[\([knt + 1]\)]*Dntmj + \ \ Dnaty[\([knt + 1]\)]\ }, 
	{j, 
        0, nd}]\)

\(natsplplot = 
      ListPlot[Table[{xnat[j], ynat[j]}, {j, 0, nd}], 
        PlotJoined\  -> \ True, 
	
        PlotRange -> {{\(-5\), 5}, {\(-4\), 4}}, 	PlotStyle -> Green, 
        AspectRatio -> 0.5, AxesLabel -> {"> x", "  ^ y"}]\)

graf[lauf] = 
      Show[natsplplot, newtonplot, ellipsplot, linienplot, punktplot, 
        Prolog  ->  AbsolutePointSize[5]]

Do[{Print["       Lauf Nummer ", l , " mit ", 
          nummer[l], " Stützpunkten " ], Show[graf[l]]}, {l, 
        1, lauf}]