SF1626 - Övning 3 - Code

Göran Andersson goeran@kth.se, 110330

Paketet SF1626 - Några bra metoder

Som bilaga till denna övning finns filen SF1626.m som innehåller kod för metoderna Taylor och CompleteSquare samt några ytterligare metoder.  Filen SF1626.nb är en minimanual för paketet. Där visas också hur du installerar paketet.

Efter att du har startat Mathematica så laddas tilläggspaketet med

Needs["SF1626`"]

Börja alltid med detta kommando om du avser att använda metoderna i paketet SF1626 . Kommandot Needs["SF1626`"] ignoreras om paketet redan är laddat.

Bakgrund

Här ges lite bakgrundsinformation för den intresserade.

Ett allmänt Taylor-polynom av grad n kan skrivas med hjälp av gradienten.

(Subscript[T, n])(𝕣)=(∑_(k=0)^n)(1/(k!))((((𝕣-(Subscript[𝕣, 0]))·∇))^k)f((Subscript[𝕣, 0]))

Metoden Taylor implementerar denna formel. Syntaxen lyder:

Taylor[ f[𝕣], 𝕣, (Subscript[𝕣, 0]), n ]

Taylor[ f[x,y,...], {x, y, ...}, {x0, y0, ...}, n ]

Exempel: Skattning runt (1,2,0) av ordning 3 för 1/(1+3(x^2)+2(y^2)+(z^2)).

Needs["SF1626`"]

ClearAll["`*"]

Taylor[(1/(1+3(x^2)+2(y^2)+(z^2))),{x,y,z},{1,2,0},3]

Expand[%]

Allmänt polynom av grad 2 för en 2D-funktion:

Taylor[f[a+h,b+k],{h,k},{0,0},2]

Här är ett annat exempel med vektornotation:

g[r_]:=((r.b)/(((1+r.a))^(3/2)))

r={x,y,z};r0={1,2,-1};
a={1,2,3};b={1,-1,2};
n=2;
pT[{x_,y_,z_}]=Expand[Taylor[g[r],r,r0,n]]

Låt oss testa noggrannheten:

Δ={0.1,-0.1,0.1};
((pT[r0+Δ]-g[r0+Δ])/(g[r0+Δ]))

Metoden CompleteSquare har följande syntax:

CompleteSquare[ p[x], x ]

metoden fungerar även för kvadratiska former. Successiv kvadratkomplettering m.a.p. x,y,... erhålls med:

CompleteSquare[ q[x, y, ...], {x, y, ...} ]

Exempel:

ClearAll["`*"]

CompleteSquare[ (x^2)+3x+5, x]

CompleteSquare[ a (t^2)+b t+c, t]

CompleteSquare[ 21 (x^2)-48 x y+43 (y^2)-30 x z+45 (z^2), {x,y,z}]

CompleteSquare[ 21 (x^2)-48 x y+43 (y^2)-30 x z+45 (z^2), {z,x,y}]

Metoden SymmetricMatrix skapar en matris som genererar polynomet q(𝕣) enligt  (Transpose[𝕣])·𝔸·𝕣.

SymmetricMatrix[ q[x, y, ...], {x, y, ...} ]

MatrixForm[
A=SymmetricMatrix[21 (x^2)-48 x y+43 (y^2)-30 x z+45 (z^2),{x,y,z}]
]

Expand[{x,y,z}.A.{x,y,z}]

Med matrisen kan olika definite-typer testas. Följande metod finns inbyggd i Mathematica.

PositiveDefiniteMatrixQ[A]

Övriga testmetoder finns i paketet SF1626, t.ex.

IndefiniteMatrixQ[SymmetricMatrix[(x^2)+3x y-(y^2),{x,y}]]

Exempel: Hitta en negativt definit symmetrisk 5×5 matris med slumpheltal.

n=5;intervall={-4,4};
rndA:=With[{m=RandomInteger[intervall,{n,n}]},Transpose[m]+m]

A=rndA;
While[!NegativeDefiniteMatrixQ[A],
A=rndA
];
MatrixForm[A]

Alla egenvärdena är negativa.

N[Eigenvalues[A]]

Några allmänna tips

En funktion kan införas på olika sätt. Exempelvis

f[x_, y_] := uttryck

innebär anrop av typ f[3,4] och

f[{x_, y_}] := uttryck

innebär anrop av typ f[{3,4}].  Vilken syntax man väljer är en smaksak.

Gradienten beräknas med

D[f, {{x1, x2, ...}}]

Hessianen beräknas med

D[f, {{x1, x2, ...},2}]

Exempel:

f[x_,y_]:=x (y^2)+(x^2)(y^3)

D[f[x,y],{{x,y}}]

D[f[x,y],{{x,y},2}]

2.28

En påminnelse: Markera en klammern till en kodcell och välj menyn Cell, Convert To, InputForm för att se koden i ren ASCII.

a)

ClearAll["`*"](* rensa egna symboler *)

f[x_,y_]:=((((x^2)+(y^2)))^n)

D[f[x,y],{{x,y}}]

2.33

ClearAll["`*"](* rensa egna symboler *)

Vi inför vektornotationer här. Double struck letters införs med escape-sekvens EscdsxEsc.

z[{x_,y_}]:=(1/100)(50-(x^2)-(y^2));
𝕩={x,y};𝕒={7,-1};𝕓={-1,3};

Gradienten blir (nabla ∇ införs med EscdelEsc).

(∇z)[{x_,y_}]=D[z[𝕩],{𝕩}]

Riktningen införs med

𝕧=((𝕓-𝕒)/(Norm[𝕓-𝕒]))

Punkten 𝕡 ligger på linjen

𝕡=𝕒+t 𝕧

Vi söker punkt där riktningsderivatan är 0.

(∇z)[𝕡].𝕧//Simplify

{t0}=t/.Solve[(∇z)[𝕡].𝕧==0]

Rasten tas därmed vid

𝕡0=𝕡/.t->t0

z[𝕡0]

Riktningsderivatan är en linjär funktion av t. Därmed antas extremvärdena vid start och slut.

{(∇z)[𝕒].𝕧,(∇z)[𝕓].𝕧}

Vi avslutar med lite grafik:

hill=Plot3D[
z[𝕩],{x,-8,8},{y,-8,8},RegionFunction->Function[{x,y,z},Norm[{x,y}]≤(Sqrt[50])],
AxesLabel->Automatic,PlotStyle->Opacity[0.7]
]

Solve[𝕡==𝕓]

t1=0;t2=4 (Sqrt[5]);

ℙ[t_]=𝕡~Join~{z[𝕡]}//Simplify

path=ParametricPlot3D[ℙ[t],{t,t1,t2},PlotStyle->Directive[Thick,Red]]

breakpoint=Graphics3D[{PointSize[0.02],Point[ℙ[t0]]}]

Show[hill,path,breakpoint]

2.37

ClearAll["`*"](* rensa egna symboler *)

T[x_,y_]:=(20/π) ArcTan[((2 Cos[x])/((E^y)-(E^(-y))))]

Gradienten blir

g[x_,y_]=D[T[x,y],{{x,y}}]//Simplify

g[(π/3),Log[2]]

Kryp i riktningen

v=((g[(π/3),Log[2]])/(Norm[g[(π/3),Log[2]]]))

Vi avslutar med lite grafik:

Plot3D[
T[x,y],{x,-(π/2),(π/2)},{y,0,3},
Mesh->{Range[0,10]},MeshFunctions->{Function[{x,y,z},z]},
ColorFunction->"TemperatureMap",AxesLabel->Automatic
]

ContourPlot[
T[x,y],{x,-(π/2),(π/2)},{y,0,3},
Contours->Range[0,10],ColorFunction->"TemperatureMap",
FrameLabel->Automatic
]

2.43

a)

ClearAll["`*"](* rensa egna symboler *)

f[{x_,y_,z_}]:=x (z^5)+x y z-9

𝕩={x,y,z};𝕩0={3,2,1};

Punkten (Subscript[𝕩, 0])=(3,2,1) ligger då på nivåytan N:f(𝕩)=0. Gradienten ∇f((Subscript[𝕩, 0])) blir:

(∇f)[{x_,y_,z_}]=D[f[𝕩],{𝕩}]

(∇f)[𝕩0]

Gradienten ∇f((Subscript[𝕩, 0])) är vinkelrät mot alla vektorer 𝕩-(Subscript[𝕩, 0]) i tangentplanet:

tangentplan=(∇f)[𝕩0].(𝕩-𝕩0)==0

Simplify[tangentplan]

b)

Vi approximerar punkten (Subscript[𝕩, 1])=(3.4,2.3,z) med punkten på tangentplanet:

z/.Solve[tangentplan/.{x->3.4,y->2.3}]

Korrekt numerisk lösning ger:

𝕩1={3.4,2.3,z};
z/.FindRoot[f[𝕩1]==0,{z,0.9}]

2.50

ClearAll["`*"]

f[{x_,y_,z_}]:=Sin[(x^2)+(y^2)]+x y z

𝕩={x,y,z};

Gradienten ∇f(𝕩) blir:

D[f[𝕩],{𝕩}]

TraditionalForm[%]

Hessianen blir:

D[f[𝕩],{𝕩,2}]

TraditionalForm[%]

2.51

ClearAll["`*"]

f[x_,y_]:=g[(x^2)-y]
de=2D[f[x,y],y]+D[f[x,y],{x,2}]+x D[f[x,y],{y,2}]==0

de2=Simplify[de/. {(x^2)-y->t}]

(g[t]/.DSolve[de2,g,t])/.t->(x^2)-y

2.53

ClearAll["`*"]

T[r_,t_]:=f[(r/(Sqrt[t]))]
de=D[T[r,t],t]==D[T[r,t],{r,2}]-(1/r)D[T[r,t],r]

de2=Simplify[de/.{r->s(Sqrt[t])},t>0]

f[s]/.DSolve[de2,f,s]

2.61

Vi använder metoden Taylor  i paketet SF1626.

Needs["SF1626`"]

a)

ClearAll["`*"]

f[x_,y_]:=(((1+x+(y^2)))^(1/2))

Utveckling av f(0+h,0+k) kring (h,k)=(0,0).

Expand[Taylor[f[0+h,0+k],{h,k},{0,0},2]]

b)

ClearAll["`*"]

f[x_,y_]:=(((x+1))^(y+1))

Utveckling av f(0+h,0+k) kring (h,k)=(0,0).

Expand[Taylor[f[0+h,0+k],{h,k},{0,0},2]]

2.62

Vi använder metoden Taylor  i paketet SF1626.

Needs["SF1626`"]

b)

ClearAll["`*"]

f[x_,y_]:=(Sqrt[1+x+y])

Utveckling av f(1+h,0+k) kring (h,k)=(0,0).

p=Expand[Taylor[f[1+h,0+k],{h,k},{0,0},2]]

För att skriva om uttrycket kan ma flytta termerna med copy and paste:

(Sqrt[2])+(h/(2 (Sqrt[2])))+(k/(2 (Sqrt[2])))+(-((h^2)/(16 (Sqrt[2])))-((h k)/(8 (Sqrt[2])))-((k^2)/(16 (Sqrt[2]))))

Använd sedan paletten  AlgebraicManipulation. Markera sedan uttrycket i sista parentesen och tryck Factor[■]. Gör därefter samma sak med term 2 och 3. Resultatet blir:

(Sqrt[2])+((h+k)/(2 (Sqrt[2])))+(-((((h+k))^2)/(16 (Sqrt[2]))))