SF1626 - Övning 5 - Code

Göran Andersson goeran@kth.se, 110407

En påminnelse: Om du vill omvandla kod till ren ASCII för att enkelt kunna mata in koden gör då följande: Markera klammern till kodcellen och välj menyn Cell, Convert To, InputForm.

Kort om syntax

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

När en metod definieras används ofta _ eller Blank i vänsterledet. Kortfattat betyder _ godtyckligt uttryck och t.ex. x_ betyder ett godtyckligt uttryck som kallas x. Hänvisning till detta i högerledet sker med x (utan _). Operatorn := betyder tilldelning av kod (fördröjd tilldelning) och = betyder tilldelning av resultat. I följande exempel tilldelas (Subscript[f, 1])(x) koden (∫_0^x)a (t^2)dx och (Subscript[f, 2])(x) resultatet (x^3). Om konstanten a ändras kommer (Subscript[f, 1]) att påverkas men inte (Subscript[f, 2]).

a=3;

f1[x_]:=Integrate[a (t^2),{t,0,x}]

f2[x_]=Integrate[a (t^2),{t,0,x}]

f1[2]

f2[2]

a=5;
{f1[2],f2[2]}

Operatorn /. eller ReplaceAll används ofta för att uföra substitutioner tillsammans med -> eller Rule. Syntaxen följer typiskt

a b+c/.{a->2,b->-3,c->4}

a b+c/.{{a->2,b->-3,c->4},{a->x,b->y,c->z}}

Lägg märke till att Solve levererar en lösning som kan användas till substitution:

Solve[(x^2)+3x+1==0]

{x1,x2}=x/.Solve[(x^2)+3x+1==0,x]

Kommandot Reduce levererar lösning i form av logiska utsagor. Lösningen kan i många fall omvandlas till den form som Solve levererar.

losning=Reduce[(x^2)+3x+1==0,x]

{ToRules[losning]}

Operatorn /@ eller Map används för att applicera funktioner på listor av data.

f/@{1,2,3}

Map[f,{{1,2,3},{4,5,6}},{2}]

Operatorn @@ eller Apply används för att applicera en funktion på en lista av argument.

arg={1,2,3};
f@@arg

Apply[f,{{1,2},{3,4},{5,6}},{1}]

Vektorer eller matriser multipliceras med . eller Dot.

{1,2,3}.{a,b,c}

{{x,y}}.{{1,2,3},{4,5,6}}

{{1,2,3},{4,5,6}}.{{x},{y},{z}}

I sistnämnda fallen ger ofta multiplikation av matris (2D) med vektor (1D) en bekvämare syntax och man behöver inte bekymra sig om transponat.

{x,y}.{{1,2,3},{4,5,6}}

{{1,2,3},{4,5,6}}.{x,y,z}

2.66

Needs["SF1626`"]

ClearAll["`*"]

Vi skall bestämma lokala extrempunkter till följande funktion.

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

Vi söker nollställen till gradienten:

𝕩={x,y};
(∇f)[{x_,y_}]=D[f[𝕩],{𝕩}]

𝕠={0,0};
{𝕩1,𝕩2}=𝕩/.Solve[(∇f)[𝕩]==𝕠,𝕩]

Den kvadratiska formen beräknas med 2:a-derivatorna (Hessianen).

ℋ[{x_,y_}]=D[f[𝕩],{𝕩,2}]

𝕙={h,k};
𝒬1[{h_,k_}]=𝕙.ℋ[𝕩1].𝕙

CompleteSquare[𝒬1[𝕙],𝕙]

Vi kan se att (Subscript[Q, 1])(𝕙) är indefinit och därmed är (Subscript[𝕩, 1]) en terrasspunkt.

𝒬2[{h_,k_}]=𝕙.ℋ[𝕩2].𝕙

CompleteSquare[𝒬2[𝕙],𝕙]

Vi kan se att (Subscript[Q, 2])(𝕙) är positivt definit och därmed är (Subscript[𝕩, 2]) en lokal maxpunkt.

Vi avslutar med grafik:

p1=𝕩1~Join~{f[𝕩1]}

p2=𝕩2~Join~{f[𝕩2]}

plot=Plot3D[
f[x,y],{x,-0.3,0.3},{y,-0.3,0.3},
PlotStyle->Opacity[0.75],MeshFunctions->{Function[{x,y,z},z]},
Mesh->{{0.001,0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.5}},
AxesLabel->Automatic,PlotPoints->75
];
points=Graphics3D[{PointSize[0.015],Blue,Point[p1],Red,Point[p2]}];
Show[plot,points]

För dig som vill veta mer

Informationen om vilken typ den kvadratiska formen är finns helt i Hessianen. Det innebär att vi inte behöver beräkna den kvadratiska formen utan kan ställa frågan “till” matrisen.

Needs["SF1626`"]

IndefiniteMatrixQ[ℋ[𝕩1]]

PositiveDefiniteMatrixQ[ℋ[𝕩2]]

4.1

ClearAll["`*"]

Mängder kan beskrivas med RegionPlot och RegionPlot3D. Här används 𝒦 (script 𝒦 som skrivs EscscKEsc) för att undvika konflikt med inbyggda K.

?K

𝒦[x_,y_]:=Abs[x]≤1&&Abs[y]<1

𝒦[.5,.3]

𝒦[1.5,.3]

RegionPlot[
𝒦[x,y],{x,-4,4},{y,-4,4},
FrameLabel->{"x","y"},
GridLines->Automatic
]

4.3

ClearAll["`*"]

Området är

𝒦[x_,y_]:=(-1≤x≤1)&&(-1≤y≤1)

Funktionen definieras:

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

Extremvärdena blir

MinValue[{f[x,y],𝒦[x,y]},{x,y}]

MaxValue[{f[x,y],𝒦[x,y]},{x,y}]

4.8

ClearAll["`*"]

Området är

𝒦[x_,y_]:=0≤y≤x≤1

Funktionen definieras:

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

Extremvärdena blir

MinValue[{f[x,y],𝒦[x,y]},{x,y}]

MaxValue[{f[x,y],𝒦[x,y]},{x,y}]

Grafik

m=Minimize[{f[x,y],𝒦[x,y]},{x,y}]

M=Maximize[{f[x,y],𝒦[x,y]},{x,y}]

p1={x,y,m[[1]]}/.m[[2]]

p2={x,y,M[[1]]}/.M[[2]]

plot1=Plot3D[
f[x,y],{x,0,2},{y,0,2},PlotRange->All,
PlotStyle->Opacity[0.75],
RegionFunction->Function[{x,y,z},𝒦[x,y]],
Mesh->{Range[0,2,0.5],Range[0,2,0.5]},
AxesLabel->{"x","y","f"}
]

plot2=Graphics3D[{
{PointSize[0.02],Red,Point[p1]},
{PointSize[0.02],Blue,Point[p2]}
}];
Show[plot1,plot2]

4.12

ClearAll["`*"]

Området och funktionen definieras:

𝒦[x_,y_]:=(x^2)+(y^2)≤4

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

RegionPlot[𝒦[x,y],{x,-2,2},{y,-2,2},
FrameLabel->Automatic,GridLines->Automatic]

Metoderna MinValue och MaxValue fugerar inte här så vi försöker manuellt. Det inre området och randen definieras:

ℐ[x_,y_]:=(x^2)+(y^2)<4

ℬ[x_,y_]:=(x^2)+(y^2)==4

Vi söker punkt där ∇f=𝟘.

innersolution={x,y}/.{ToRules[Reduce[{(∂_({{x,y}}))f[x,y]=={0,0},ℐ[x,y]}]]}

Andra möjligheter finns på randen.

g[y_]=Simplify[f[x,y],ℬ[x,y]]

Reduce[g'[y]==0]

Detta ger två punkter, {(-2,0),(2,0)}. Sammanlagt har vi följande punkter att undersöka:

points={{-2,0},{2,0}}~Join~innersolution

Funktionens värde i punkterna blir:

values=Apply[f,points,{1}]

Därmed kan det minsta och det största värdet beräknas.

{m,M}={Min[values],Max[values]}

Värdena antas i punkterna (0,0) samt (0,±1).

Grafik

Vi kan få numerisk lösning:

m=NMinimize[{f[x,y],𝒦[x,y]},{x,y}]//Chop

M=NMaximize[{f[x,y],𝒦[x,y]},{x,y}]//Chop

Funktionen är jämn m.a.p. både x och y så om (x,y)=(0,1) är en lösning är även (0,-1) det.

p1={x,y,m[[1]]}/.m[[2]]

p2={x,y,M[[1]]}/.M[[2]]

p3={x,-y,M[[1]]}/.M[[2]]

plot=Plot3D[
f[x,y],{x,-2,2},{y,-2,2},PlotRange->All,
PlotStyle->Opacity[0.75],
RegionFunction->Function[{x,y,z},𝒦[x,y]],
Mesh->{Range[-2,2],Range[-2,2]},
AxesLabel->Automatic
];
points=Graphics3D[{
{PointSize[0.02],Red,Point[p1]},
{PointSize[0.02],Blue,Point[p2],Point[p3]}
}];
Show[plot,points]

4.25

ClearAll["`*"]

Bivillkoret är

g[x_,y_]:=13(x^2)+13(y^2)+10x y-72

plot1=ContourPlot[g[x,y]==0,{x,-3,3},{y,-3,3},FrameLabel->{"x","y"},GridLines->Automatic]

Funktionen definieras:

d[x_,y_]:=(Sqrt[(x^2)+(y^2)])

Extremvärdena blir

MinValue[{d[x,y],g[x,y]==0},{x,y}]

MaxValue[{d[x,y],g[x,y]==0},{x,y}]

Manuell lösning

Alt. 1

solution=Reduce[{
(∂_({{x,y}}))d[x,y]==λ (∂_({{x,y}}))g[x,y],
g[x,y]==0
},{x,y}]

{x,y,d[x,y]}/.{ToRules[solution]}

Alt. 2

G={(∂_({{x,y}}))d[x,y],(∂_({{x,y}}))g[x,y]}

solution=Reduce[{
Det[G]==0,
g[x,y]==0
},{x,y}]

Grafik

m=Minimize[{d[x,y],g[x,y]==0},{x,y}]

M=Maximize[{d[x,y],g[x,y]==0},{x,y}]

Metoderna returnerar endast en punkt där globalt extremvärde antas. I detta fall kan vi på grund av symmetri sluta oss till att även ±(x,y) är möjliga lösningar. Dvs. värdena antas i punkterna ±(Sqrt[2])(1,1) samt ±(3/(Sqrt[2]))(1,-1).

p1={x,y,m[[1]]}/.m[[2]]

p2={-x,-y,m[[1]]}/.m[[2]]

p3={x,y,M[[1]]}/.M[[2]]

p4={-x,-y,M[[1]]}/.M[[2]]

plot1=Plot3D[
d[x,y],{x,-3,3},{y,-3,3},PlotRange->All,
PlotStyle->Opacity[0.75],Mesh->{Range[-3,3],Range[-3,3]},
AxesLabel->{"x","y","d"}
];
plot2=ContourPlot3D[
g[x,y]==0,{x,-3,5},{y,-3,3},{z,0,3},
ContourStyle->Opacity[0.25],
Mesh->None
];
plot3=Graphics3D[{
{PointSize[0.02],Red,Point[p1],Point[p2]},
{PointSize[0.02],Green,Point[p3],Point[p4]}
}];
Show[plot1,plot2,plot3]

Parametrisering:

sol=Solve[g[r Cos[t],r Sin[t]]==0,r]//Simplify

R=r/.sol[[2]]

𝕣[t_]={R Cos[t],R Sin[t],R}

plot4=ParametricPlot3D[
𝕣[t],{t,0,2π},
PlotRange->{All,All,{0,4}},
AxesLabel->{"x","y","d"}
];
Show[plot4,plot3]

4.29

ClearAll["`*"]

Bivillkoret är g(x,y)=0

g[x_,y_,z_]:=((((x/a)))^2)+((((y/b)))^2)+((((z/c)))^2)-1

Funktionen definieras:

V[x_,y_,z_]:=8x y z

positive:=(x>0&&y>0&&z>0&&a>0&&b>0&&c>0)

Maximum blir

MaxValue[{V[x,y,z],g[x,y,z]==0&&positive},{x,y,z}]

4.32

ClearAll["`*"]

Bivillkoren är g(x,y,z)=0∧h(x,y,z)=0

g[x_,y_,z_]:=(x^2)+(y^2)+(z^2)-2

h[x_,y_,z_]:=(x^2)+(y^2)-z

plot1=ContourPlot3D[{g[x,y,z]==0,h[x,y,z]==0},{x,-2,2},{y,-2,2},{z,-2,2},AxesLabel->{"x","y","z"},Mesh->False,ContourStyle->{Opacity[0.5],Opacity[0.5]},PlotPoints->25]

Funktionen definieras:

f[x_,y_,z_]:=x+y+z

Maxvärdet blir:

MaxValue[{f[x,y,z],g[x,y,z]==0&&h[x,y,z]==0},{x,y,z}]

Grafik

M=Maximize[{f[x,y,z],g[x,y,z]==0&&h[x,y,z]==0},{x,y,z}]

p0={x,y,z}/.M[[2]]

plot2=ContourPlot3D[
f[x,y,z],{x,-5,5},{y,-5,5},{z,-5,5},
Mesh->None,Contours->{M[[1]]//N},
ContourStyle->{
Directive[Red,Opacity[0.5]],
Directive[Green,Opacity[0.5]]
}
];
plot3=Graphics3D[{PointSize[0.015],Point[p0]}];
Show[plot1,plot2,plot3]