<Text-field style="Heading 1" layout="Heading 1">Init</Text-field>

Bedste rette linje / bedste plan

Rev. 29.08.25, stema@dtu.dk

restart; with(linalg): with(plots):with(plottools):

with(LinearAlgebra):

with(VectorCalculus):

BasisFormat(false):

with(combinat):

base2d:= proc(x,y, len) description "2d basis plot"; ########################################################### local e1_2d,mark1_2d, e2_2d, mark2_2d, axlong1, axlong2, Nlen; ########################################################### if nargs=2 then Nlen:= [-1.5,1.5, -1.5,1.5] else Nlen:= len end if; e1_2d:= arrow([0,0], [1,0], 0.05, 0.10, 0.15, color=red): e2_2d:= arrow([0,0], [0,1], 0.05, 0.10, 0.15, color=blue): mark1_2d:= x-> textplot([1.2, -0.2,`x`], color=black, font=[TIMES, ITALIC, 20]): mark2_2d:= y-> textplot([-0.2, 1.2,`y`], color=black, font=[TIMES, ITALIC, 20]): axlong1:= plot([u,0, u=Nlen[1]...Nlen[2]], color=black); axlong2:= plot([0,u, u=Nlen[3]...Nlen[4]], color=black); display([e1_2d, mark1_2d(x), e2_2d, mark2_2d(y), axlong1, axlong2], scaling=constrained, axes=none): ########################################################### end: ###########################################################

<Text-field style="Heading 1" layout="Heading 1">Et s\303\246t punkter i planen</Text-field>

Antal punkter: N

N:= 7;

pt:= [[-.2783858614, .1781784561], [.3079459675, 1.038621939], [.6959791960, 1.250528672], [-1.254507455, .1688706500], [-.7752164941, -.1132856455], [1.106758651, .8110377852], [1.011027796, 1.696848490]];

ptset:= seq(disk(pt[i], 0.03, color=red), i=1...N):

display(ptset, base2d(x,y));

<Text-field style="Heading 1" layout="Heading 1">En bev\303\246gelig test-linje</Text-field>

r:= (u,v,x0, y0)-> [x0+u*cos(v), y0+ u*sin(v)];

Afstand (med fortegn) fra (x, y) til linjen med punkt (x0, y0) p\303\245 linjen og vinkel v:

f:= (x,y,v,c)-> sin(v)*x-cos(v)*y + c;

solve(f(x0+u*cos(v), y0+ u*sin(v), v, c)=0, c);

Afstanden udtrykt via vilk\303\245rligt punkt (x0, y0) p\303\245 linjen:

f:= (x,y,v,x0, y0)-> simplify(sin(v)*x-cos(v)*y -sin(v)*x0 + cos(v)*y0);

Afstands-cirkelskiver:

disking:= (v,x0, y0)-> seq(disk(pt[i], evalf(f(op(pt[i]), v, x0, y0)), color=gray), i=1...7):

linBev:= (v,x0,y0)-> display(base2d(x,y), ptset, plot([op(r(u, v, x0,y0)), u=-3...3], color=blue, thickness=2));

Med eksplicit valg (g\303\246t) af (x0, y0) og v f\303\245s:

display(ptset, disking(Pi/3, 1/2,1), linBev(Pi/3, 1/2,1),scaling=constrained);

Generelt med variabel (x0, y0) og v:

Sview:= (v,x0,y0)-> display(ptset, disking(v, x0,y0), linBev(v, x0,y0));

Forskellige valg af (x0, y0) og v giver displays:

Sview(Pi/3, 1/2,1);

Sview(Pi/8, 1/2,1);

Sview(Pi, 1/2,1);

<Text-field style="Heading 1" layout="Heading 1">Massemidtpunkt Cm og rotation om Cm</Text-field>

Cm:= (1/N)*add(pt[i], i=1...N);

CmPt:= disk(Cm, 0.07, color=green):

Den ovenfor g\303\246ttede linje g\303\245r ikke igennem det gr\303\270nne massemidtpunkt:

display(CmPt, Sview(Pi/3, 1/2,1));

Benyt massemidtpunkt som nyt linjepunkt:

display(CmPt,Sview(Pi/3, op(Cm)));

<Text-field style="Heading 1" layout="Heading 1">Den direkte analyse med henblik p\303\245 at finde den bedste vinkel</Text-field>

v:= 'v';

A, B, C og Inertimomentmatricen:

A:= add((pt[i][2]-Cm[2])^2, i=1...N);

B:= add((pt[i][1]-Cm[1])^2, i=1...N);

C:= add((pt[i][1]-Cm[1])*(pt[i][2]-Cm[2]), i=1...N);

Imom:= Matrix(2,2, [[A, -C], [-C, B]]);

Fejlen:

S:=v-> A*cos(v)^2 + B*sin(v)^2 - 2*C*cos(v)*sin(v); 'S(v)' = S(v);

Check: S(v) er summen af kvadratafstandene:

simplify(add(f(op(pt[i]), v,x0,y0)^2, i=1..N), expand)

Differentialkvotientne af S(v):

SP:= unapply(diff(S(v), v),v);

Bedste vinkel:

v0:= fsolve(SP(v)=0, v=0...evalf(Pi/2));

Check op mod arctan udtrykket i noten:

evalf(0.5*arctan(2*C/(B-A)));

SVD alternativet Extract til projektioner

<Text-field style="Heading 1" layout="Heading 1">7 point example </Text-field>

N:= nops(pt):

Cm:= (1/N)*add(pt[i], i=1...N):

M:= Transpose(Matrix([seq(pt[i]-Cm, i=1...N)]));

Sigma:= SingularValues(M, output = 'list');

SigmaMat:= (eps1, eps2)-> DiagonalMatrix([Sigma[1]*eps1, Sigma[2]*eps2], 2, N);

SigmaMat(1,1);

U, Vt := SingularValues(M, output = ['U', 'Vt']):

projV:= convert(Vt[1], list);

e1:= convert(Transpose(U)[1], list);

ghostProj:= [seq(Cm + Sigma[1]*projV[i]*e1, i=1...N)];

ghostProjFig:= display(seq(circle(ghostProj[i], 0.05, color=green, thickness=3), i=1...N)):

theta:= argument(U[1,1]+ I*U[2,1]);

SVD projektionsoperator:

projSVD:= (i, eps1, eps2)-> convert(evalm(Cm+Transpose(U.SigmaMat(eps1, eps2).Vt)[i]), list);

projSVDFig:= (eps1, eps2)-> seq(circle(projSVD(i, eps1,eps2), 0.1, color=gray), i=1...N):

SVD projektioner:

ptset:= seq(disk(pt[i], 0.03, color=red), i=1...N):

CmPt:= disk(Cm, 0.07, color=green):

disking:= (v,x0, y0)-> seq(disk(pt[i], evalf(f(op(pt[i]), v, x0, y0)), color=gray), i=1...N):

linBev:= (v,x0,y0)-> display(ptset, plot([op(r(u, v, x0,y0)), u=-3...3], color=blue, thickness=4)):

Sview:= (v,x0,y0)-> display(ptset, disking(v, x0,y0), linBev(v, x0,y0)):

SviewPure:= (v,x0,y0)-> display(ptset, linBev(v, x0,y0)):

SviewPurePure:= (v,x0,y0)-> plot([op(r(u, v, x0,y0)), u=-3...3], color=blue, thickness=4):

Den bedste linje og den v\303\246rste:

projTo1:= t-> display(ptset, CmPt, SviewPure(theta, op(Cm)), projSVDFig(1,t), gridlines, labels=["x", "y"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], view=[-3...4, -1...4], scaling=constrained);

projTo2:= t-> display(ptset, CmPt, SviewPure(theta+Pi/2, op(Cm)), projSVDFig(t,1), gridlines, labels=["x", "y"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], view=[-3...4, -1...4], scaling=constrained);

display(projTo1(0), ghostProjFig);

display(seq(projTo1(1-q/10), q=0..10), insequence=true);

display(seq(projTo2(1-q/10), q=0..10), insequence=true);

Direkte ortogonale projektioner:

for i from 1 to N do proj1[i]:= Cm + [cos(theta), sin(theta)]*DotProduct(Vector(pt[i]-Cm), Vector([cos(theta), sin(theta)])); proj1Pt[i]:= circle(proj1[i], 0.06, color=red, thickness=3); proj2[i]:= Cm + [cos(theta+Pi/2), sin(theta+Pi/2)]*DotProduct(Vector(pt[i]-Cm), Vector([cos(theta+Pi/2), sin(theta+Pi/2)])); proj2Pt[i]:= circle(proj2[i], 0.06, color=red, thickness=3); end do:

projOrto1:= display(SviewPure(theta, op(Cm)), seq(proj1Pt[i], i=1...N), scaling=constrained);

projOrto2:= display(SviewPure(theta+Pi/2, op(Cm)), seq(proj2Pt[i], i=1...N), scaling=constrained);

SVD projektionen er samme som ortogonale projektioner:

seq(projSVD(i, 1, 0)-proj1[i], i=1...N);

<Text-field style="Heading 1" layout="Heading 1">Eksempel 4; Udvidet med Projektions-overvejelser</Text-field>

rol:= evalf(rand(-5...5)/10):

pt:= evalf([seq([1+i/20 + rol(), 1.5+2*(-(i/20))], i=-10...10)]):

N:= nops(pt):

Cm:= (1/N)*add(pt[i], i=1...N):

M:= Transpose(Matrix([seq(pt[i]-Cm, i=1...N)]));

Sigma:= SingularValues(M, output = 'list');

SigmaMat:= (eps1, eps2)-> DiagonalMatrix([Sigma[1]*eps1, Sigma[2]*eps2], 2, N);

SigmaMat(1,1);

U, Vt := SingularValues(M, output = ['U', 'Vt']):

theta:= argument(U[1,1]+ I*U[2,1]);

SVD projektionsoperator:

projSVD:= (i, eps1, eps2)-> convert(evalm(Cm+Transpose(U.SigmaMat(eps1, eps2).Vt)[i]), list);

projSVDFig:= (eps1, eps2)-> seq(circle(projSVD(i, eps1,eps2), 0.1, color=gray), i=1...N):

SVD projektioner:

ptset:= seq(disk(pt[i], 0.03, color=red), i=1...N):

CmPt:= disk(Cm, 0.07, color=green):

disking:= (v,x0, y0)-> seq(disk(pt[i], evalf(f(op(pt[i]), v, x0, y0)), color=gray), i=1...N):

linBev:= (v,x0,y0)-> display(ptset, plot([op(r(u, v, x0,y0)), u=-3...3], color=blue, thickness=4)):

Sview:= (v,x0,y0)-> display(ptset, disking(v, x0,y0), linBev(v, x0,y0)):

SviewPure:= (v,x0,y0)-> display(ptset, linBev(v, x0,y0)):

SviewPurePure:= (v,x0,y0)-> plot([op(r(u, v, x0,y0)), u=-3...3], color=blue, thickness=4):

Den bedste linje og den v\303\246rste:

display(seq(projTo1(1-q/10), q=0..10), insequence=true);

display(seq(projTo2(1-q/10), q=0..10), insequence=true);

Direkte ortogonale projektioner:

projOrto1:= display(SviewPure(theta, op(Cm)), seq(proj1Pt[i], i=1...N), scaling=constrained);

projOrto2:= display(SviewPure(theta+Pi/2, op(Cm)), seq(proj2Pt[i], i=1...N), scaling=constrained);

SVD projektionen er samme som ortogonale projektioner:

seq(projSVD(i, 1, 0)-proj1[i], i=1...N);

<Text-field style="Heading 1" layout="Heading 1">Det tredimensionale tilf\303\246lde: Projektionernes identitet</Text-field>

rol:= evalf(rand(-5...5)/20):

pt:= evalf([seq([1+i/20 + rol(), 1.5+2*(-(i*rol()/10)), 3*i/20 + rol()], i=-10...10)]):

N:= nops(pt):

ptFig:= display(seq(sphere(pt[i], 0.06, color=red, style=patchnogrid), i=1...N), scaling=constrained);

Fig:= display(ptFig, axesfont=["HELVETICA", 20], labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

viewing:= [-1...3, 0.5...2.5, -4...4 ]

display(ptFig,view=viewing);

Cm:= (1/N)*add(pt[i], i=1...N):

CmPt:= sphere(Cm, 0.1, color=green, style=patchnogrid):

M:= Transpose(Matrix([seq(pt[i]-Cm, i=1...N)]));

Sigma:= SingularValues(M, output = 'list'):

U, Vt := SingularValues(M, output = ['U','Vt']):

e1:= convert(Transpose(U)[1], list);

e2:= convert(Transpose(U)[2], list);

e3:= convert(Transpose(U)[3], list);

Sigma;

SigmaMat:= (eps1, eps2, eps3)-> DiagonalMatrix([Sigma[1]*eps1, Sigma[2]*eps2, Sigma[3]*eps3], 3, N);

SigmaMat(1,1,1);

U, Vt := SingularValues(M, output = ['U', 'Vt']):

ax1:= t-> convert(evalm(Cm+t*e1), list);ax1(1);

ax1Fig:= spacecurve(ax1(t), t=-4...4, color=blue, thickness=4):

ax2:= t-> convert(evalm(Cm+t*e2), list);ax2(1);

ax2Fig:= spacecurve(ax2(t), t=-2...2, color=khaki, thickness=4):

ax3:= t-> convert(evalm(Cm+t*e3), list);ax1(1);

ax3Fig:= spacecurve(ax3(t), t=-2...2, color=black, thickness=4):

Fig:= display(ptFig, CmPt,ax1Fig, ax2Fig, ax3Fig, axesfont=["HELVETICA", 20], labelfont = ["HELVETICA", 20], tickmarks=[4,4,4], view=viewing);

display(Fig);

SVD projektionsoperator

projRumSVD:= (i, eps1, eps2, eps3)-> convert(evalm(Cm+Transpose(U.SigmaMat(eps1, eps2, eps3).Vt)[i]), list);fr

projRumSVDFig:= (eps1, eps2, eps3)-> display(seq(sphere(projRumSVD(i, eps1,eps2, eps3), 0.06, color=gray, style=patchnogrid), i=1...N), scaling=constrained):

projRumSVDFig(1,0,0);

SVD projektioner til principale akser:

projRumTo1:= t-> display(ptFig, CmPt, ax1Fig, projRumSVDFig(1,t,t),labels=["x", "y", "z"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], scaling=constrained);

projRumTo1(0.);

projRumTo2:= t-> display(ptFig, CmPt, ax2Fig, projRumSVDFig(t,1,t),labels=["x", "y", "z"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], scaling=constrained);

projRumTo2(0.);

projRumTo3:= t-> display(ptFig, CmPt, ax3Fig, projRumSVDFig(t,t,1),labels=["x", "y", "z"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], scaling=constrained);

projRumTo3(0.5);

SVD projektioner til principale planer:

plan23:= plot3d(evalm(Cm +u*e2 + v*e3), u=-2..2, v=-0.9...0.9, color=cyan, grid=[3,3]):

plan13:= plot3d(evalm(Cm +u*e1 + v*e3), u=-2..2, v=-0.9...0.9, color=cyan, grid=[3,3]):

plan12:= plot3d(evalm(Cm +u*e1 + v*e2), u=-2..2, v=-0.9...0.9, color=cyan, grid=[3,3]):

projRumTo23:= t-> display(ptFig, CmPt, ax2Fig, ax3Fig, plan23, projRumSVDFig(t,1,1),labels=["x", "y", "z"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], scaling=constrained);

projRumTo23(0.);

projRumTo13:= t-> display(ptFig, CmPt, ax1Fig, ax3Fig, plan13, projRumSVDFig(1,t,1),labels=["x", "y", "z"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], scaling=constrained);

projRumTo13(0.);

projRumTo12:= t-> display(ptFig, CmPt, ax1Fig, ax2Fig, plan12, projRumSVDFig(1,1,t),labels=["x", "y", "z"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], scaling=constrained);

projRumTo12(2.0);

Projektionerne direkte ortogonalt

for i from 1 to N do proj1[i]:= Cm + e1*DotProduct(Vector(pt[i]-Cm), Vector(e1)); proj1Pt[i]:= sphere(proj1[i], 0.06, color=red, thickness=3, style=patchnogrid); proj2[i]:= Cm + e2*DotProduct(Vector(pt[i]-Cm), Vector(e2)); proj2Pt[i]:= sphere(proj2[i], 0.06, color=red, thickness=3, style=patchnogrid); proj3[i]:= Cm + e3*DotProduct(Vector(pt[i]-Cm), Vector(e3)); proj3Pt[i]:= sphere(proj3[i], 0.06, color=red, thickness=3,style=patchnogrid); end do:

Fig:= display(seq(proj1Pt[i], i=1...N), scaling=constrained,ax1Fig, CmPt, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+21, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

Fig:= display(seq(proj2Pt[i], i=1...N), scaling=constrained,ax2Fig,CmPt, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+22, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

Fig:= display(seq(proj3Pt[i], i=1...N), scaling=constrained,ax3Fig,CmPt, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+23, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

De tre planprojektioner og tilh\303\270rende fejl:

From scratch fra plan-normal:

pm11:= add((pt[i][1]-Cm[1])^2, i=1...N);

pm22:= add((pt[i][2]-Cm[2])^2, i=1...N);

pm33:= add((pt[i][3]-Cm[3])^2, i=1...N);

pm12:= add((pt[i][1]-Cm[1])*(pt[i][2]-Cm[2]), i=1...N);

pm13:= add((pt[i][1]-Cm[1])*(pt[i][3]-Cm[3]), i=1...N);

pm23:= add((pt[i][2]-Cm[2])*(pt[i][3]-Cm[3]), i=1...N);

pm21:= pm12;

pm31:= pm13;

pm32:= pm23;

Pmom:= Matrix(3,3, [[pm11, pm12, pm13], [pm21, pm22,pm23], [pm31, pm32, pm33]]);

Eigenvectors(Pmom);

Sigma[1]^2;Sigma[2]^2;Sigma[3]^2;

e1; e2; e3;

for i from 1 to N do projP23[i]:= Cm + (pt[i]-Cm)- e1*DotProduct(Vector(pt[i]-Cm), Vector(e1)); projP23Pt[i]:= sphere(projP23[i], 0.06, color=red, thickness=3, style=patchnogrid); projP13[i]:= Cm + (pt[i]-Cm)-e2*DotProduct(Vector(pt[i]-Cm), Vector(e2)); projP13Pt[i]:= sphere(projP13[i], 0.06, color=red, thickness=3, style=patchnogrid); projP12[i]:= Cm +(pt[i]-Cm)- e3*DotProduct(Vector(pt[i]-Cm), Vector(e3)); projP12Pt[i]:= sphere(projP12[i], 0.06, color=red, thickness=3,style=patchnogrid); end do:

plan23:= plot3d(evalm(Cm +u*e3 + v*e2), u=-2..2, v=-0.9...0.9, color=cyan, grid=[3,3]):

projP23Fig:= display(seq(projP23Pt[i], i=1...N), scaling=constrained);

Fig:= display(projP23Fig, plan23,ax3Fig, ax2Fig, ax1Fig, CmPt, view=viewing, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+24, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

plan13:= plot3d(evalm(Cm +u*e3 + v*e1), u=-1...1, v=-3...3, color=cyan, grid=[3,3]):

projP13Fig:= display(seq(projP13Pt[i], i=1...N), scaling=constrained);

Fig:= display(projP13Fig, plan13, ax3Fig, ax2Fig, ax1Fig, CmPt, view=viewing, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+25, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

plan12:= plot3d(evalm(Cm +u*e2 + v*e1), u=-0.9...0.9, v=-3...3, color=cyan, grid=[3,3]):

projP12Fig:= display(seq(projP12Pt[i], i=1...N), scaling=constrained);

Fig:= display(projP12Fig, plan12,ax3Fig, ax2Fig, ax1Fig, CmPt, view=viewing, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+26, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

Fig:= display(ptFig, plan13, ax3Fig, ax2Fig, ax1Fig, CmPt, view=viewing, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+27, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; #interface(plotdevice=default);

Fig:= display(ptFig, plan23, ax3Fig, ax2Fig, ax1Fig, CmPt, view=viewing, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+28, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; #interface(plotdevice=default);

Fig:= display(ptFig, plan12, ax3Fig, ax2Fig, ax1Fig, CmPt, view=viewing, axesfont=["HELVETICA", 20], view=viewing, labelfont = ["HELVETICA", 20], tickmarks=[4,4,4]);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVD", convert(100+29, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; #interface(plotdevice=default);

Klynger

<Text-field style="Heading 1" layout="Heading 1">Deling af de 7 punkter i to disjunkte punkts\303\246t:</Text-field>

Brute force partitioning:

#pt:= [[-2,-1], [-2, 1], [0, -1], [0, 1], [2, -1], [2,1]];

#rol:= evalf(rand(-20...20)/20):

#pt:= [seq([rol(), rol()], a=1...10)];

N:= nops(pt);

Cm:= (1/N)*add(pt[i], i=1...N):

FM0:= add(Norm(Vector(pt[j]- Cm))^2, j=1....N);

S:= [seq(i, i=1...N)];

par:= setpartition(S): parN:= nops(par);

seq(nops(par[i]), i=1...parN):

q:= 1: for i from 1 to parN do if nops(par[i])=2 and nops(par[i][1])>1 and nops(par[i][2])>1 then par2[q]:= par[i]; q:=q+1; else end if; end do:

par2set:= seq(par2[j], j=1...q-1):

P1:= i-> [seq(pt[par2set[i][1][w]],w=1...nops(par2set[i][1]))]:

P1Fig:= i-> seq(disk(P1(i)[j], 0.05, color=red), j=1...nops(par2set[i][1]));

P2:= j-> [seq(pt[par2set[j][2][w]],w=1...nops(par2set[j][2]))]:

P2Fig:= i-> seq(disk(P2(i)[j], 0.05, color=blue), j=1...nops(par2set[i][2]));

P1(2);

mu1:= i-> (1/nops(P1(i)))*add(P1(i)[j], j=1...nops(P1(i))): mu1(3);

m1Fig:= i-> circle(mu1(i), 0.08, color=red):

mu2:= i-> (1/nops(P2(i)))*add(P2(i)[j], j=1...nops(P2(i))): mu2(3);

m2Fig:= i-> circle(mu2(i), 0.08, color=blue):

FM1:= i-> add(Norm(Vector(P1(i)[j]- mu1(i)))^2, j=1....nops(P1(i)));

FM2:= i-> add(Norm(Vector(P2(i)[j]- mu2(i)))^2, j=1....nops(P2(i)));

sort([seq([FM1(i) + FM2(i), i], i=1...q-1)])[1...5];

minim:= sort([seq([FM1(i) + FM2(i), i], i=1...q-1)])[1][2];

FM1(minim);

FM2(minim);

ptset:= seq(disk(pt[i], 0.02, color=gray), i=1...N):

figRaw:=display(ptset, P1Fig(minim), P2Fig(minim),gridlines, labels=["x", "y"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], view=[-2...2, -1...2], scaling=constrained);

interface(plotdevice=jpeg, plotoutput=cat("SVDklynge", convert(100+1, string), ".jpeg"), plotoptions = "height=2000,width=2000"); FigRaw; #interface(plotdevice=default);

Linje-separation p\303\245 f\303\270rste principale akse (repr\303\246senteret p\303\245 x-aksen med nummerering):

Cm:= (1/N)*add(pt[i], i=1...N):

M:= Transpose(Matrix([seq(pt[i]-Cm, i=1...N)]));

Sigma:= SingularValues(M, output = 'list');

SigmaMat:= (eps1, eps2)-> DiagonalMatrix([Sigma[1]*eps1, Sigma[2]*eps2], 2, N);

SigmaMat(1,1);

U, Vt := SingularValues(M, output = ['U', 'Vt']):

projV:= convert(Vt[1], list);

cand:= sort([seq([Sigma[1]*projV[i], i], i=1...N)]);

Kandidat-sekvens:

seq(cand[q][2], q=1....N);

Delinger:

for q from 2 to N-2 do P1new[q]:= [seq(pt[cand[w][2]], w=1...q)]: P2new[q]:= [seq(pt[cand[u][2]], u=q+1...N)]: end do:

mu1new:= i->(1/nops(P1new[i]))*add(P1new[i][j], j=1...nops(P1new[i])): mu1new(3);

mu2new:= i-> (1/nops(P2new[i]))*add(P2new[i][j], j=1...nops(P2new[i])): mu2new(3);

FM1new:= i-> add(Norm(Vector(P1new[i][j]- mu1new(i)))^2, j=1....nops(P1new[i])): FM1new(2);

FM2new:= i-> add(Norm(Vector(P2new[i][j]- mu2new(i)))^2, j=1....nops(P2new[i])): FM2new(2);

sort([seq([FM1new(i) + FM2new(i), i], i=1...q-1)])[1...5];

minim:= sort([seq([FM1new(i) + FM2new(i), i], i=1...q-1)])[1][2];

ptset:= seq(disk(pt[i], 0.02, color=gray), i=1...N):

P1newFig:= i-> seq(circle(P1new[i][j], 0.07, color=red), j=1...nops(P1new[i]));

P2newFig:= i-> seq(circle(P2new[i][j], 0.07, color=blue), j=1...nops(P2new[i]));

display( figRaw, gridlines, labels=["x", "y"],axesfont=["HELVETICA", 20],labelfont = ["HELVETICA", 20], view=[-2...2, -1...2], scaling=constrained);

Fig:= display(P1newFig(3), P2newFig(3), figRaw, scaling=constrained, gridlines, labels=["x", "y"],axesfont=["HELVETICA", 40],labelfont = ["HELVETICA", 40], view=[-2...2, -1...2], scaling=constrained);

interface(plotdevice=jpeg, plotoutput=cat("SVDklynge", convert(100+2, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

<Text-field style="Heading 1" layout="Heading 1">Deling af random stor punktm\303\246ngde med SVD projektion i 2D</Text-field>

rol:= evalf(rand(-20...20)/20):

pt:= [seq([rol(), rol()], a=1...45)]:

N:= nops(pt);

Cm:= (1/N)*add(pt[i], i=1...N):

Linje-separation p\303\245 f\303\270rste principale akse (repr\303\246senteret p\303\245 x-aksen med nummerering):

Cm:= (1/N)*add(pt[i], i=1...N):

M:= Transpose(Matrix([seq(pt[i]-Cm, i=1...N)]));

Sigma:= SingularValues(M, output = 'list');

SigmaMat:= (eps1, eps2)-> DiagonalMatrix([Sigma[1]*eps1, Sigma[2]*eps2], 2, N);

SigmaMat(1,1);

U, Vt := SingularValues(M, output = ['U', 'Vt']):

projV:= convert(Vt[1], list):

cand:= sort([seq([Sigma[1]*projV[i], i], i=1...N)]):

Kandidat-sekvens:

seq(cand[q][2], q=1....N);

Delinger:

for q from 2 to N-2 do P1new[q]:= [seq(pt[cand[w][2]], w=1...q)]: P2new[q]:= [seq(pt[cand[u][2]], u=q+1...N)]: end do:

mu1new:= i->(1/nops(P1new[i]))*add(P1new[i][j], j=1...nops(P1new[i])): mu1new(3);

mu2new:= i-> (1/nops(P2new[i]))*add(P2new[i][j], j=1...nops(P2new[i])): mu2new(3);

FM1new:= i-> add(Norm(Vector(P1new[i][j]- mu1new(i)))^2, j=1....nops(P1new[i])): FM1new(2);

FM2new:= i-> add(Norm(Vector(P2new[i][j]- mu2new(i)))^2, j=1....nops(P2new[i])): FM2new(2);

sort([seq([FM1new(i) + FM2new(i), i], i=2...N-2)])[1...4];

minimNew:= sort([seq([FM1new(i) + FM2new(i), i], i=2...N-2)])[1][2];

ptset:= seq(disk(pt[i], 0.02, color=gray), i=1...N):

P1newFig:= i-> seq(circle(P1new[i][j], 0.03, color=red, thickness=4), j=1...nops(P1new[i]));

P2newFig:= i-> seq(circle(P2new[i][j], 0.03, color=blue, thickness=4), j=1...nops(P2new[i]));

klyngeCm1:= disk(mu1new(minimNew), 0.04, color=gray):

klyngeCm2:= disk(mu2new(minimNew), 0.04, color=gray):

Fig:= display(P1newFig(minimNew), P2newFig(minimNew), klyngeCm1, klyngeCm2, scaling=constrained, gridlines, labels=["x", "y"],axesfont=["HELVETICA", 40],labelfont = ["HELVETICA", 40], view=[-1.2...1.2, -1.2...1.2], scaling=constrained);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVDklynge", convert(100+2, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

<Text-field style="Heading 1" layout="Heading 1">Deling af random stor punktm\303\246ngde med SVD projektion i 3D</Text-field>

rol:= evalf(rand(-20...20)/20):

pt:= [seq([rol(), rol(), rol()], a=1...45)]:

N:= nops(pt);

Cm:= (1/N)*add(pt[i], i=1...N);

Linje-separation p\303\245 f\303\270rste principale akse (repr\303\246senteret p\303\245 x-aksen med nummerering):

M:= Transpose(Matrix([seq(pt[i]-Cm, i=1...N)]));

Sigma:= SingularValues(M, output = 'list');

SigmaMat:= (eps1, eps2)-> DiagonalMatrix([Sigma[1]*eps1, Sigma[2]*eps2], 2, N);

SigmaMat(1,1);

U, Vt := SingularValues(M, output = ['U', 'Vt']):

projV:= convert(Vt[1], list):

cand:= sort([seq([Sigma[1]*projV[i], i], i=1...N)]):

Kandidat-sekvens:

seq(cand[q][2], q=1....N);

Delinger:

for q from 2 to N-2 do P1new[q]:= [seq(pt[cand[w][2]], w=1...q)]: P2new[q]:= [seq(pt[cand[u][2]], u=q+1...N)]: end do:

mu1new:= i->(1/nops(P1new[i]))*add(P1new[i][j], j=1...nops(P1new[i])): mu1new(3);

mu2new:= i-> (1/nops(P2new[i]))*add(P2new[i][j], j=1...nops(P2new[i])): mu2new(3);

FM1new:= i-> add(Norm(Vector(P1new[i][j]- mu1new(i)))^2, j=1....nops(P1new[i])): FM1new(2);

FM2new:= i-> add(Norm(Vector(P2new[i][j]- mu2new(i)))^2, j=1....nops(P2new[i])): FM2new(2);

sort([seq([FM1new(i) + FM2new(i), i], i=2...N-2)])[1...4];

minimNew:= sort([seq([FM1new(i) + FM2new(i), i], i=2...N-2)])[1][2];

ptset:= seq(sphere(pt[i], 0.02, color=gray, style=patchnogrid), i=1...N):

P1newFig:= i-> seq(sphere(P1new[i][j], 0.04, color=red, style=patchnogrid), j=1...nops(P1new[i]));

P2newFig:= i-> seq(sphere(P2new[i][j], 0.04, color=blue, style=patchnogrid), j=1...nops(P2new[i]));

klyngeCm1:= sphere(mu1new(minimNew), 0.07, color=gray,style=patchnogrid):

klyngeCm2:= sphere(mu2new(minimNew), 0.07, color=gray, style=patchnogrid):

Fig:= display(P1newFig(minimNew), P2newFig(minimNew), klyngeCm1, klyngeCm2, labels=["x", "y", "z"],axesfont=["HELVETICA", 40],labelfont = ["HELVETICA", 40], scaling=constrained);

display(Fig);

#interface(plotdevice=jpeg, plotoutput=cat("SVDklynge", convert(100+4, string), ".jpeg"), plotoptions = "height=2000,width=2000"); Fig; interface(plotdevice=default);

LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbW9HRiQ2LlEifkYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0Y0LyUqc3ltbWV0cmljR0Y0LyUobGFyZ2VvcEdGNC8lLm1vdmFibGVsaW1pdHNHRjQvJSdhY2NlbnRHRjQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZDLyUpZXF1YWxpdHlHRjQvJStleGVjdXRhYmxlR0Y0Ri8=