7  Løsbare ligningsystem

I svært mange sammenhenger trenger man å løse et lineært ligningssystem.
Det betyr at vi har flere ukjente variabler som ganges med konstanter og summeres.

Et generelt eksempel er

$$ \[\begin{align} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n &= b_1\\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n &= b_2 \\ \vdots \qquad & \\ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n &= b_n \end{align}\] $$

Dette kan skrives kompakt på matriseform:

\[ A x = b \]

I dette kapittelet fokuserer vi først på det kvadratiske tilfellet \(n \times n\), der systemet lar seg løse entydig, som vil være tilfellet om, \(\det A \neq 0\), som også vil si at det finnes en inversmatrise \(A^{-1}\) .

7.0.1 Løse lineærsystemet \(A x = b\) med Julia

La \[ A = \begin{bmatrix} 2 & 1 & 0 & 0 \\ 1 & 3 & 1 & 0 \\ 0 & 1 & 2 & 1 \\ 0 & 0 & 1 & 2 \end{bmatrix}, \qquad b = \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix}. \] Med bruk av Julia vil vi

  1. Løse for \(x\) ved å bruke inversen \(A^{-1}b\) .
  2. Finne løsningen direkte med å bruke \(A \backslash b\) .
  3. Verifisere at \(A^{-1}A = I\) og \(AA^{-1} = I\).
# Definer A og b
A = [2 1 0 0;
     1 3 1 0;
     0 1 2 1;
     0 0 1 2]

b = [1, 2, 3, 4]

# 1) Løsning via invers
Ainv = inv(A)
x_inv = Ainv * b

# 2) Direkte løsning (anbefalt i praksis)
x_dir = A \ b

# Skriv ut løsninger og verifiser likhet
println("x via inv(A)*b =\n", x_inv)
println("x via A\\b     =\n", x_dir)

# 3) Verifisering av identitet
I_left  = Ainv * A
I_right = A * Ainv

println("\nAinv * A =\n", I_left)
println("\nA * Ainv =\n", I_right)

Om man kun skal løse en enkelt ligning \(A x = b\), så finner man sjelden selve inversen. Det er mer effektivt å løse \(A x = b\) direkte, som I Julia gøres med back-slash, A\b, som kan “tolkes” som A under brøkstreken. Merk at for “matematiske tekseter” bør notasjonen \(A \backslash b\) generelt sett unngås.

7.1 Anvendelser

Vi vil nå ta for oss eksempler på slike ligningssystem som oppstår i teknologi.

7.1.1 Balansere kjemiske reaksjoner

Ligninger av typen \(Ax=b\) oppstår når vi skal prøve å å løse kjemiske problemer.

7.1.1.1 Forbrenning av glukose

Vi ser på (ubalansert) forbrenning av glukose:

\[ C_6H_{12}O_6 + O_2 \;\;\to\;\; CO_2 + H_2O. \]

Vi ønsker å finne koeffisientene \(a, b, c, d\) slik at

\[ a\,C_6H_{12}O_6 + b\,O_2 \;\;\to\;\; c\,CO_2 + d\,H_2O \]

blir en balansert kjemisk reaksjon. Dette betyr at antall karbon-, hydrogen- og oksygenatomer er det samme på begge sider.

Antall atomer i hvert stoff

Vi setter opp en tabell der radene tilsvarer grunnstoffer og kolonnene tilsvarer forbindelser.
Hvert innslag er antall atomer av det aktuelle grunnstoffet i ett molekyl.

Grunnstoff $C_6H_{12}O_6$ $O_2$ $CO_2$ $H_2O$
C 6 0 1 0
H 12 0 0 2
O 6 2 2 1

Lingingsystem

Reaktanter får positive kolonner, produkter får negative kolonner.
Dette gir balansematrisen

\[ M = \begin{bmatrix} 6 & 0 & -1 & 0 \\ 12 & 0 & 0 & -2 \\ 6 & 2 & -2 & -1 \end{bmatrix}. \] La

\[ x = \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} \]

være de ukjente koeffisientene, balanserer vi problemet med å løse

\[ M x = 0. \]

Systemet har uendelig mange løsninger, siden vi bare bestemmer forholdet mellom \(a,b,c,d\). Derfor setter vi \(a=1\).

Dette svarer til

\[ M \begin{bmatrix} 1 \\ b \\ c \\ d \end{bmatrix} = 0. \]

Når vi setter \(a=1\), kan vi skrive systemet slik:

\[ M \begin{bmatrix} 1 \\ b \\ c \\ d \end{bmatrix} = 1 \cdot m_1 + b \cdot m_2 + c \cdot m_3 + d \cdot m_4 = 0, \]

der \(m_j\) er kolonnene i matrisen \(M\):

\[ m_1 = \begin{bmatrix} 6 \\ 12 \\ 6 \end{bmatrix}, \quad m_2 = \begin{bmatrix} 0 \\ 0 \\ 2 \end{bmatrix}, \quad m_3 = \begin{bmatrix} -1 \\ 0 \\ -2 \end{bmatrix}, \quad m_4 = \begin{bmatrix} 0 \\ -2 \\ -1 \end{bmatrix}. \]

Dette kan da uttrykkes som

\[ b\,m_2 + c\,m_3 + d\,m_4 = -m_1. \] eller

\[ \begin{bmatrix} m_2 & m_3 & m_4 \end{bmatrix} x' = -m_1 \]

hvor \(x' = (b, c, d)\)

using LinearAlgebra

# Kolonnene m2, m3, m4
M = [ 0  -1   0;
      0   0  -2;
      2  -2  -1 ]

# Høyresiden -m1
mm1 = [-6, -12, -6]

# Løs ligningssystemet
x = M \ mm1   # (b, c, d)
println("b, c, d = ", x)

Løsningen er \(b=6,\; c=6,\; d=6\), eller

\[ C_6H_{12}O_6 + 6\,O_2 \;\;\to\;\; 6\,CO_2 + 6\,H_2O. \]

7.1.1.2 Vanadiumredoksbatteri

Batterier basert på vandium-redoks-reaksjoner er en lovende teknologi for storskala energilagring.
Her må vi balansere både atomer og ladning.

Den totale celle-reaksjonen er

\[ V^{2+} + VO_2^+ + H^+ \;\;\to\;\; V^{3+} + VO^{2+} + H_2O. \]

Vi skriver dette mer generelt med koeffisienter

\[ a\,V^{2+} + b\,VO_2^+ + c\,H^+ \;\;\to\;\; d\,V^{3+} + e\,VO^{2+} + f\,H_2O. \]

Antall atomer og ladning i hvert stoff

Størrelse \(V^{2+}\) \(VO_2^+\) \(H^+\) \(V^{3+}\) \(VO^{2+}\) \(H_2O\)
V 1 1 0 1 1 0
O 0 2 0 0 1 1
H 0 0 1 0 0 2
Ladning +2 +1 +1 +3 +2 0

Balansematrise

\[ M = \begin{bmatrix} 1 & 1 & 0 & -1 & -1 & 0 \\ 0 & 2 & 0 & 0 & -1 & -1 \\ 0 & 0 & 1 & 0 & 0 & -2 \\ 2 & 1 & 1 & -3 & -2 & 0 \end{bmatrix}, \qquad x = \begin{bmatrix} a \\ b \\ c \\ d \\ e \\ f \end{bmatrix}. \]

Balansen betyr \(Mx=0\).

Feste én variabel

Velg \(a=1\). Da kan vi skrive

\[ M \begin{bmatrix} 1 \\ b \\ c \\ d \\ e \\ f \end{bmatrix} = 0 \;\;\;\Longrightarrow\;\;\; M_{:,2:}\,x' = -m_1, \]

der \(x' = (b,c,d,e,f)\) og \(m_1\) er første kolonne i \(M\).

A = [1 0 -1 -1 0;
     2 0  0 -1 -1;
     0 1  0  0 -2;
     1 1 -3 -2  0]   # M[:, 2:end]

b = -[1, 0, 0, 2]    # -m1

x = A \ b
println("b, c, d, e, f = ", x)

Løsningen her er ikke i heltall, men i \(1/3\), vi kan enten skalere, eller så kan vi, f.eks. sette \(a=3\), som svarer til

A = [1 0 -1 -1 0;
     2 0  0 -1 -1;
     0 1  0  0 -2;
     1 1 -3 -2  0]   # M[:, 2:end]

b = -3*[1, 0, 0, 2]    # -m1

x = A \ b
println("b, c, d, e, f = ", x)

Dette gir:

\[ 3\,V^{2+} + 2\,VO_2^+ + 6\,H^+ \;\;\to\;\; 4\,V^{3+} + \,VO^{2+} + 3\,H_2O. \]

7.2 Interpolasjon

Mange former for maskinlæring, spesielt veiledet læring (supervised learning) kan forstås som en interpolasjon. Dette gjelder spesielt om dataene ikke er støyete, som f.eks. om simuleringer brukes som treningsdata. Mange former for interpolasjon krever at man løser ikke-lineære problemer, men her skal vi ta for oss interpolasjon som kan løses lineært, altså det involverer å løse problemet å løse problemet \(Ax = b\).

En stor fordel med slik interpolasjon (utover at det er raskt) er at man er garantert en løsning når \(A\) er invertibel, og at dette er gunstig når operasjonen skal automatiseres eller brukes som et grunnlag for mer avanserte beregninger.

7.2.0.1 Polynominterpolasjon

En vanlig måte å representere funksjoner er ved hjelp av polymer, f.eks.

\[ p(x) = c_0 + c_1 x + c_2 x^2 + c_3 x^3 + c_4 x^4 + \ldots \]

osv. Denne formen dukker f.eks. opp i en Taylor ekspansjon. Slike funksjoner er enkle å derivere og integrere, og selve funksjons kan uttrykkes som en indreprodukt (her til 4 orden)

\[ p(x) = \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \\ c_4 \end{bmatrix}^\top \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \\ x^4 \end{bmatrix} \]

Polynominterpolasjon har den fordelen at det kan uttrykkes som et lineært interpolasjonsproblem. Vi ønsker å interpolere funksjonen

\[ f(x) = x e^{-x} \]

basert på 4 målepunkter og kunnskapen om at \(f(0)=0\), så vi kan sette konstantleddet lik null på forhånd. Om vi lar målepunktene være \(x_1, x_2, x_3, x_4\) og krever at \(p(x_i) = f(x_i)\) for alle \(i\)., så gir dette problemet

\[ \begin{bmatrix} x_1 & x_1^2 & x_1^3 & x_1^4 \\ x_2 & x_2^2 & x_2^3 & x_2^4 \\ x_3 & x_3^2 & x_3^3 & x_3^4 \\ x_4 & x_4^2 & x_4^3 & x_4^4 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ c_4 \end{bmatrix} = \begin{bmatrix} f(x_1) \\ f(x_2) \\ f(x_3) \\ f(x_4) \end{bmatrix}. \]

Matrisen på denne formen kalles en Vandermonde-matrise

For eksempel, med punktene \(x = 1, 2, 3, 4\), får vi

\[ \begin{bmatrix} 1 & 1 & 1 & 1 \\ 2 & 4 & 8 & 16 \\ 3 & 9 & 27 & 81 \\ 4 & 16 & 64 & 256 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ c_4 \end{bmatrix} = \begin{bmatrix} f(1) \\ f(2) \\ f(3) \\ f(4) \end{bmatrix}. \]

Julia eksempel:

using Plots

f(x) = x*exp(-x)
xs = 1:4
b = f.(xs)

# Vandermonde-matrise uten konstantledd
A = [x^j for x in xs, j in 1:4]

# Løs systemet
c = A \ b
println("c = ", c)

# Definer interpolasjonspolynomet
p(x) = c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4

# Plot funksjonen og interpolasjonen
plot(f, 0, 4, label="f(x) = x e^{-x}", lw=2)
plot!(p, 0, 4, label="interpolasjonspolynom", lw=2, color="red")
scatter!(xs, b, label="datapunkter", ms=6, color="red")

Resultatet er som følger:

7.2.1 Padé-interpolasjon

Polynominterpolasjon kan ofte gi dårlig oppførsel utenfor interpolasjons-området, f.eks. divergerer det alltid for store verdier av \(x\).

Et ofte oversett alternativ er å bruke en Padé-approksimasjon, som kan uttrykkes rasjonell funksjon

\[ r(x) = \frac{p(x)}{q(x)} = \frac{c_0 + c_1 x + \cdots + c_m x^m}{1 + d_1 x + \cdots + d_n x^n}. \]

Ved å velge polynomgrader \(m\) og \(n\) får vi totalt \(m+n+1\) ukjente koeffisienter.
Disse kan bestemmes ved å kreve at \(r(x_i) = f(x_i)\) for et tilsvarende antall punkter \(x_1, \dots, x_{m+n+1}\).

Dette gir ligninger på formen

\[ c_0 + c_1 x_i + \cdots + c_m x_i^m = f(x_i)\,\big(1 + d_1 x_i + \cdots + d_n x_i^n\big), \quad i = 1, \dots, m+n+1. \]

Som et lineært system kan dette skrives

\[ \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^m & -f(x_1)x_1 & -f(x_1)x_1^2 & \cdots & -f(x_1)x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^m & -f(x_2)x_2 & -f(x_2)x_2^2 & \cdots & -f(x_2)x_2^n \\ \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{m+n+1} & x_{m+n+1}^2 & \cdots & x_{m+n+1}^m & -f(x_{m+n+1})x_{m+n+1} & -f(x_{m+n+1})x_{m+n+1}^2 & \cdots & -f(x_{m+n+1})x_{m+n+1}^n \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_m \\ d_1 \\ d_2 \\ \vdots \\ d_n \end{bmatrix} = \begin{bmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_{m+n+1}) \end{bmatrix}. \]

Dermed kan Padé-approksimasjon løses ved å sette opp et \(Ax=b\)-problem, helt analogt med polynominterpolasjon, men med både teller- og nevnerkoeffisienter som ukjente.

For mange former for funksjoner kan vi kjenne til oppførselen eller formen for små og store \(x\). For vår test-funksjonen

\[ f(x) = x e^{-x}, \]

så kan vi anta at vi vet at \(f(0) = 0\) og at \(f(x) \to 0\) når \(x \to \infty\).

Vi kan da. f.eks. velge en Padé-form som både går gjennom origo og har innebygd avtagende oppførsel:

\[ r(x) = \frac{c_1 x}{1 + d_1 x + d_2 x^2}. \]

  • Teller starter med \(x\), så \(r(0)=0\).
  • Nevneren vokser som \(x^2\) for store \(x\), slik at \(r(x) \to 0\).

Vi har da ligningene

\[ c_1 x_i = f(x_i)\,\big(1 + d_1 x_i + d_2 x_i^2\big), \quad i=1,2,3. \]

Eller på matriseform:

\[ \begin{bmatrix} x_1 & -f(x_1)x_1 & -f(x_1)x_1^2 \\ x_2 & -f(x_2)x_2 & -f(x_2)x_2^2 \\ x_3 & -f(x_3)x_3 & -f(x_3)x_3^2 \end{bmatrix} \begin{bmatrix} c_1 \\ d_1 \\ d_2 \end{bmatrix} = \begin{bmatrix} f(x_1) \\ f(x_2) \\ f(x_3) \end{bmatrix}. \]

Julia-eksempel:

using Plots

f(x) = x*exp(-x)
xs = [0.5, 1, 3]
b = f.(xs)

# Bygg opp Padé-systemet
A = [[x  -f(x)*x  -f(x)*x^2] for x in xs]
A = vcat(A...)

u = A \ b
c1, d1, d2 = u
println("c1=$c1, d1=$d1, d2=$d2")

r(x) = (c1*x) / (1 + d1*x + d2*x^2)

plot(f, 0, 5, label="f(x) = x e^{-x}", lw=2)
plot!(r, 0, 5, label="Padé-approksimasjon", lw=2, color="red", ls=:dash)
scatter!(xs, b, label="målepunkter", ms=6, color="red")

Her har vi valgt 3 interpolasjonspunkter, bedre fit kan man kanskje få med å velge 3 andre. Men dette er uansett ikke spesielt imponerende, så vi prøver til høyere order

\[ r(x) = \frac{c_1 x + c_2 x^2}{1 + d_1 x + d_2 x^2 + d_3 x^3}. \]

using Plots

f(x) = x*exp(-x)

# Choose 5 interpolation points
xs = [0.5, 1, 2, 3, 4]
b = f.(xs)

# Build Padé system (numerator: c1*x + c2*x^2, denominator: 1 + d1*x + d2*x^2 + d3*x^3)
A = [[x  x^2  -f(x)*x  -f(x)*x^2  -f(x)*x^3] for x in xs] 
A = vcat(A...)   # splat the rows into vcat

u = A \ b 
c1, c2, d1, d2, d3 = u 
println("c1=$c1, c2=$c2, d1=$d1, d2=$d2, d3=$d3")

r(x) = (c1*x + c2*x^2) / (1 + d1*x + d2*x^2 + d3*x^3)

plot(f, 0, 5, label="f(x) = x e^{-x}", lw=2)
plot!(r, 0, 5, label="Padé-approksimasjon", lw=2, color="red", ls=:dash)
scatter!(xs, b, label="målepunkter", ms=6, color="red")

Dette gir

Dette er gir et svært godt fit.

7.3 Interpolasjon med basisfunksjoner

Et generelt prinsipp er å bygge en funksjon som en lineær kombinasjon av basisfunksjoner:

\[ f(x) \approx a_1 \varphi_1(x) + a_2 \varphi_2(x) + \cdots + a_m \varphi_m(x). \]

Her er \(\varphi_j(x)\) valgt på forhånd, og de ukjente er koeffisientene \(a_j\).

7.3.1 Eksempel: Gaussiske basisfunksjoner

Velg to basisfunksjoner

\[ \varphi_1(x) = e^{-(x-0)^2}, \qquad \varphi_2(x) = e^{-(x-1)^2}. \]

Vi ønsker at \(f(x_i) = y_i\) i to punkter \((x_1,y_1)\) og \((x_2,y_2)\).

Dette gir et lineært system

\[ \begin{bmatrix} \varphi_1(x_1) & \varphi_2(x_1) \\ \varphi_1(x_2) & \varphi_2(x_2) \end{bmatrix} \begin{bmatrix} a_1 \\ a_2 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix}. \]

7.4 Differensial-ligninger

7.4.1 Poisson’s ligninger i 1D

I kapptellet om dynamiske systemer, så vi på fordelinger som utvikler seg over tid.
Ofte er vi i stedet interessert i stasjonære tilstander – løsninger som oppfyller en balanse mellom en indre respons og en ekstern kilde.

Et viktig eksempel er Poisson-ligningen i 1D:

\[ - \frac{d^2 u}{dx^2} = f(x). \]

Her er \(f(x)\) en kjent kildefunksjon, mens \(u(x)\) er løsningen vi vil finne.

Situasjoner som gir denne typen ligning: - Varmeledning: \(u(x)=T(x)\) er temperatur, \(q(x)/\kappa\) er varmeproduksjon, hvor
- Elektrostatikk: \(u(x)=\phi(x)\) er potensial, \(f(x)=\rho(x)/\varepsilon\) er ladningstetthet.
- Mekanikk: \(u(x)\) er forskyvning, \(f(x)\) er ytre last.
- Diffusjon: \(u(x)=c(x)\) er konsentrasjon, \(f(x)=s(x)/D\) er kilde- eller slukrate.

For å løse ligningen trenger vi grensebetingelser. En vanlig situasjon er Dirichlet-betingelser:

\[ u(0)=u(1)=0. \]

For eksempel kan en oppvarmet stang ha endene holdt ved konstant temperatur.

Vi deler intervallet inn i \(N\) punkter og approksimerer den andrederiverte med endelige differanser. Da får vi et lineært ligningssystem

\[ A u_{\mathrm{vec}} = f_{\mathrm{vec}}\,, \]

der \(A\) er en tridiagonal Laplacian matrise (som tidligere). For eksempel, med \(N=5\) indre punkter:

\[ A = \frac{1}{h^2} \begin{bmatrix} 2 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 2 \end{bmatrix}, \]

der \(h=1/(N+1)\) er punktavstanden.

7.4.1.1 Eksempel: Varmeledning

Vi kan ta for seg et problem som ligner på det som vi jobbet med i forrige kapittel, med en stang som er holdt ved konstant temperatur ved sidene.

Men denne ganger vi interessert i en jevn varmekilde på midten (istedet for en initialbetingelse) og vi er interessert i temperaturfordelingen etter systemet har termalisert, det vil si at tidsavhengigheten har forsvunnet og vi får en konstant temperaturfordeling.

Som tidligere jobber vi med “effektive” enheter, hvor temperaturen ved kanten på 0 kan forstås i forhold til en konstant temperatur (som kan være 0 grader Celsisus), lengden på platen er satt til 1, slik at koordinater er relativt til den totale lengden.

Om vi setter \(f(x)=q(x)/\kappa\), og diskretiserer problemet kan vi løse dette direkte som et ligningssystem

using LinearAlgebra
using Plots

function main()
    # --- Parameters ---
    L   = 1.0           # Lengde på stangen
    n   = 100            # Antall indre punkter (ukjente)
    κ   = 1.0            # Varmeledningsevne
    dx  = L / (n + 1)    # Romlig steg

    # --- Grid (indre punkter) ---
    xgrid = (1:n) .* dx

    # --- Systemmatrix (Dirichlet: u=0 i endene) ---
    main_diag = fill( 2.0, n)
    off_diag  = fill(-1.0, n-1)
    A = (κ/dx^2) * Tridiagonal(off_diag, main_diag, off_diag)

    # --- Varmeproduksjon q(x) ---
    # f.eks. en jevn kilde i midten av stangen
    q = zeros(n)
    for i in 1:n
        xi = xgrid[i]
        if 0.4 <= xi <= 0.6
            q[i] = 1.0
        end
    end

    # --- Høyreside (Poisson: -T'' = q/κ -> A*T = q) ---
    f = q

    # --- Løs systemet A*u = f ---
    T = A \ f

    # --- Plot resultat ---
    plot(
        xgrid, T,
        lw = 2,
        xlabel = "x",
        ylabel = "Temperatur T(x)",
        title  = "Stasjonær temperaturfordeling med indre varmekilde",
        label  = false,
        color  = :red
    )
    savefig("poisson_heat.png")
    display(current())
end

main()

Om du utforsker ulike verdier av \(\kappa\), \(q\), eller endrer lengden på stangen, vil du se at kurvene får samme form. Løsningen av Poisson-ligningen har en innebygd skalainvarians: dersom både varmekilden \(q\) og varmeledningsevnen \(\kappa\) multipliseres med samme faktor, endres ikke formen på løsningen – bare temperaturnivået. Slike strukturelle egenskaper ligger til grunn for mange empiriske ingeniørlov­er. Ofte er slike lover funnet ut fra eksperiment, men de kan like gjerne komme fra simuleringer. Effektive kurverepresentasjoner kan man for eksempel finne ved “smart” interpolasjon, som vi har sett tidligere.

7.4.1.2 Eksempel 2: Hengebro

Vi skal nå se på et annet problem fra mekanikk.
Vi ser for oss en tynn og stram vaier – for eksempel en liten hengebro eller line for linedans– som er festet i begge ender og bærer sin egen vekt. I tillegg har vi to personer som går på broen/linen.

I dette tilfellet vil buen på broen bestemmes av en balanse mellom spenning i vaieren og ytre krefter (tyngdekraft og personenes vekt).
Hvis vi lar \(u(x)\) beskrive nedbøyningen av broen, får vi ligningen

\[ -\,T\,\frac{d^2 u}{dx^2} = p(x), \]

der
- \(T\) er strekkspenningen (konstant langs vaieren),
- \(p(x)\) er den påførte kraften per lengdeenhet (tyngdekraft + personer).

Randbetingelsene er at vaieren er festet i endene:

\[ u(0) = u(1) = 0. \]

Tyngdekraften virker som en jevn nedoverlast \(p_g(x) = -p_0\), mens hver person representeres som en lokalisert last (smalt intervall i \(x\)):

\[ p(x) = -p_0 - F_1(x) - F_2(x), \]

der

\[ F_i(x) = \begin{cases} f_i, & x_i - \tfrac{w}{2} \le x \le x_i + \tfrac{w}{2},\\[4pt] 0, & \text{ellers.} \end{cases} \]

Her er \(f_i\) størrelsen på kraften (vekt per lengdeenhet), \(x_i\) posisjonen til personen, og \(w\) bredden på kontaktområdet.

Etter diskretisering av andrederivert med endelige differanser får vi som før et lineært ligningssystem

\[ A\,u = p, \]

der \(A\) er en tridiagonal matrise som representerer den diskrete Laplace-operatoren, og \(p\) er en vektor som inneholder summen av tyngdekraft og personlastene. Når vi løser dette systemet, får vi nedbøyningen \(u(x)\) langs hele vaieren.

Julia-kode:

using LinearAlgebra, Plots

# --- Grid and parameters ---
L  = 1.0 
n  = 400 
h  = L / (n + 1)
x  = (1:n) .* h

# --- Matrix:  -T u'' = p  ->  A*u = p ---
Tt = 1.0 
main_diag = fill(2.0, n)
off_diag  = fill(-1.0, n - 1)
A = (Tt / h^2) * Tridiagonal(off_diag, main_diag, off_diag)

# --- Loads ---
p = -0.2 .* ones(n)   # gravity (uniform downward load)

w  = 0.05              # footprint width
x1, F1 = 0.10, 2.8     # person 1 position and force
x2, F2 = 0.68, 2.6     # person 2 position and force

p .-= F1 .* ((x1 - w/2 .<= x) .& (x .<= x1 + w/2))
p .-= F2 .* ((x2 - w/2 .<= x) .& (x .<= x2 + w/2))

# --- Solve and plot ---
u = A \ p 

plot(x, u, lw=2, label="u(x): deflection",
     xlabel="x", ylabel="u",
     title="Suspension-like wire: gravity + two walkers")
plot!(x, 0.05 .* p, lw=1, ls=:dash, label="scaled load (downward < 0)")

Enkel oppgave: Simuler hva som ville skjedd om gravitasjonen virket oppover istedet. Dette kan ses på som en rekke mindre snorer som holdt hengebroen oppe ovenfra, som i større hengebroer.

7.4.1.3 Eksempel 3: Tau trukket over ujevnt underlag.

Tenk deg et tau som er strukket mellom to punkter over et ujevnt, men mykt underlag.
Vi ønsker at tauet skal følge overflaten så tett som mulig.
Hvis det oppstår små gliper mellom tauet og underlaget, kan vi tenke oss at det er festet små fjærer som trekker tauet svakt ned mot overflaten.
Fjærene gir en kraft som er proporsjonal med avstanden mellom tau og underlag, altså omtrent som \(\lambda\,(u(x) - f(x))\),
der \(f(x)\) beskriver formen på underlaget.

Spenningen i tauet motvirker krumning, mens fjærene forsøker å trekke tauet mot overflaten.
Dette gir ligningen

\[ -\,T\,\frac{d^2 u}{dx^2} + \lambda\,(u(x) - f(x)) = 0, \]

der \(T\) er strekkspenningen i tauet, og \(\lambda\) beskriver hvor stivt eller ettergivende underlaget er.

En slik modell kalles en Winkler-overflate eller et elastisk fundament.
Hver punkt på underlaget modelleres som en uavhengig fjær, som presser tilbake når noe legger seg oppå. Dette er en idealisert modell som kan brukes (Eksemplene er fra ChatGPT og kan være usikre):

  • for å beskrive jernbanesviller og spor som ligger på elastiske masser,
  • for rørledninger og kabler som hviler på eller under jord,
  • for bjelker og plater som ligger på isolerende eller dempende lag,
  • og til og med for proteser, madrasser og skumkompositter der kontakt med et mykt underlag må modelleres.

Mer informasjon om Winkler metoden kan f.eks. i denne artikkelen: https://pubs.rsc.org/en/content/articlelanding/2018/sm/c7sm02062g

Etter diskretisering får vi et lineært system der den kontinuerlige

\[ -\,T\,\frac{d^2 u}{dx^2} + \lambda\,(u - f) = 0 \]

erstattes av \[ \bigl(\lambda I + T A \bigr)u = \lambda f, \]

der er en tridiagonal Laplace matrise

\[ A = \frac{1}{h^2} \begin{bmatrix} 2 & -1 & 0 & \cdots \\ -1 & 2 & -1 & \cdots \\ 0 & -1 & 2 & \ddots \\ \vdots & \vdots & \ddots & \ddots \end{bmatrix}. \]

Ved å dele begge sider på \(T\) kan vi skrive systemet i en dimensjonsløs form

\[ (A + kI)\,u = kf, \]

der \(k = \lambda / T\) angir hvor sterkt tauet presses mot underlaget.
For store \(k\) følger tauet overflaten tett; for små \(k\) dominerer strekkspenningen, og tauet blir jevnere.

Vi vil nå ta for oss en ru overflate som vi beskriver med tilfeldige tall. Selve overflaten beskriver vi her med en enkel parabolisk form

using Random, LinearAlgebra, Plots

# --- Grid ---
L = 1.0
n = 400
h = L / (n + 1)
x = (1:n) .* h


Random.seed!(123)

# --- Parabolic "surface" + noise ---
a = -2.0
b = 0.5
c = 0.5

f_base = a .* (x .- b).^2 .+ c

noise_level = 0.03

f = f_base .+ noise_level .* randn(n)

# --- Rope over soft surface: (k I + A) u = k f ---
# k controls how tightly the rope follows the surface
k = 700

main_diag = fill(2.0, n)
off_diag  = fill(-1.0, n - 1)
A = (1 / h^2) * Tridiagonal(off_diag, main_diag, off_diag)

M = k * I + A
u = M \ (k .* f)

# --- Plot ---
plot(x, f, lw=1, label="ru overflate f(x)")
plot!(x, f_base, lw=2, ls=:dash, label="glatt parabel f_base(x)")
plot!(x, u, lw=2, label="tau-form u(x)",
     xlabel="x", ylabel="høyde",

7.4.1.4 Eksempel 4: Utjevning av støyete data

I eksempelet med tauet/vaieren over et mykt ruglete underlag, jevner tauet ut små ujevnheter i underlaget, men følger de store strukturene. I dataanalyse kan vi bruke nøyaktig samme prinsipp for å fjerne støy fra måledata: vi ønsker en glatt kurve som følger hovedformen, men ignorerer små tilfeldige variasjoner.

Vi starter med den samme funksjonen som vi har sett på tidligere

\[ f(x) = x e^{-x}, \]

og legger til tilfeldig støy for å simulere målefeil.
Deretter bruker vi en utjevningsmodell som ligner på tauet i strekk:

\[ (\lambda I + A)\,u = \lambda f_{\text{noisy}}, \]

der \(A\) er Laplace-operatoren (samme som tidligere), og \(\lambda\) bestemmer hvor sterkt vi ønsker at kurven \(u\) skal følge de støyete dataene.

  • Stor \(\lambda\): følger data tett, lite glatting
  • Liten \(\lambda\): jevner kraftig ut, mister detaljer

Løsningen \(u\) blir en balanse mellom glatthet og tilpasning, på samme måte som tauet fant en balanse mellom strekk og underlag.

using LinearAlgebra, Random, Plots

# Funksjon og støy
n = 200 
xs = range(0, 5, length=n)
f_true = xs .* exp.(-xs)

Random.seed!(1)
f_noisy = f_true .+ 0.02*randn(n)

# Bygg 1D Laplace-operator (indre punkter, med "naturlige" randbetingelser)
h = step(xs)
main_diag = fill(2.0, n)
off_diag  = fill(-1.0, n-1)
A = Tridiagonal(off_diag, main_diag, off_diag) / h^2 

# "Screened Poisson" smoothing: (λI + A) u = λ f_noisy
λ = 100.0   # tuning parameter
M = λ*I + A 
u_smooth = M \ (λ * f_noisy)

# Plot resultat
plot(xs, f_noisy, alpha=0.4, label="Støyete data", lw=1, color=:red)
plot!(xs, u_smooth, lw=2, label="Poisson-glattet (λ=$λ)")
plot!(xs, f_true, lw=2, ls=:dash, color=:black, label="Sann funksjon")
xlabel!("x"); ylabel!("f(x)")
title!("Poisson-basert utjevning av støyete data")

Poisson-utjevning (eller screened Poisson smoothing) gir en enkel og robust måte å glatte data på. I stedet for å bruke statistiske filtre, finner vi en glatt funksjon ved å løse et lineært ligningssystem. Den samme ligningen brukes i mange felt. Innen dataanalyse og maskinlæring kalles metoden Tikhonov-regulering og er tett knyttet til ridge regresjon.

Legg merke til at selv om vi her bare løser et lineært system av typen \(A u = b\),
kan ligningen også sees som resultatet av å minimere et kompromiss mellom to størrelser –
hvor jevnheten i \(u\) balanseres mot hvor tett den følger \(f\).
Dette ligger veldig tett på hva vi skal lære om i minste-kvadraters metode.

7.5 Løsninger av \(A x = b\)

Vi har nå sett flere eksempler på lineære systemer i anvendelser, men la oss se nærmere på hvordan ulike typer matriser kan løses effektivt.
Forskjellige matrisestrukturer gir helt ulike metoder og geometriske tolkninger.

7.5.1 Diagonal matrise

Dersom \(D\) er en diagonal matrise, \(D = \mathrm{diag}(d)\), der \(d\) er en vektor med diagonale elementer \(d_i\),
får vi et svært enkelt system:

\[ A x = b \quad \Rightarrow \quad d_i x_i = b_i. \]

Løsningen finnes elementvis:

\[ x_i = \frac{b_i}{d_i}. \]

For en diagonal matrise er altså inversen også enkel:
\(D^{-1} = \mathrm{Diag}(1/d)\).

7.5.2 Triangulær matrise

En triangulær matrise har nuller enten over eller under diagonale, som ofte skrives som \(U\) og \(L\).

Om vi ser på systemet

\[ U x = b, \qquad U = \begin{bmatrix} u_{11} & u_{12} & u_{13}\\ 0 & u_{22} & u_{23}\\ 0 & 0 & u_{33} \end{bmatrix}. \]

Skriver vi ligningene ut komponentvis får vi

\[ \begin{eqnarray} u_{11}x_1 + u_{12}x_2 + u_{13}x_3 = b_1,\\ u_{22}x_2 + u_{23}x_3 = b_2,\\ u_{33}x_3 = b_3. \end{eqnarray} \]

Her kan vi starte nederst: den siste ligningen inneholder bare \(x_3\), som vi dermed kan finne direkte:

\[ x_3 = \frac{b_3}{u_{33}}. \]

Setter vi dette inn i den nest siste ligningen, kan vi løse \(x_2\), og deretter \(x_1\):

\[ x_2 = \frac{b_2 - u_{23}x_3}{u_{22}}, \qquad x_1 = \frac{b_1 - u_{12}x_2 - u_{13}x_3}{u_{11}}. \]

Dette kalles tilbakesubstitusjon (back substitution).
Algoritmen er enkel og effektiv, og krever omtrent \(n^2/2\) multiplikasjoner for en \(n\times n\) matrise.

7.5.2.1 Nedre triangulær matrise

For en nedre triangulær matrise

\[ L = \begin{bmatrix} \ell_{11} & 0 & 0\\ \ell_{21} & \ell_{22} & 0\\ \ell_{31} & \ell_{32} & \ell_{33} \end{bmatrix}, \]

får vi systemet

\[ \begin{eqnarray} \ell_{11}x_1 = b_1,\\ \ell_{21}x_1 + \ell_{22}x_2 = b_2,\\ \ell_{31}x_1 + \ell_{32}x_2 + \ell_{33}x_3 = b_3. \end{eqnarray} \]

Her starter vi ovenfra – først finner vi \(x_1\), deretter \(x_2\), osv.:

\[ x_1 = \frac{b_1}{\ell_{11}}, \qquad x_2 = \frac{b_2 - \ell_{21}x_1}{\ell_{22}}, \qquad x_3 = \frac{b_3 - \ell_{31}x_1 - \ell_{32}x_2}{\ell_{33}}. \]

Dette kalles fremoversubstitusjon (forward substitution).

7.5.2.2 Julia-eksempel

using LinearAlgebra

U = [3.0  1.0  -1.0;
     0.0  2.0   1.0;
     0.0  0.0   4.0]

b = [2.0, 1.0, 8.0]

x = UpperTriangular(U) \ b

Å løse et generelt lineært system \(A x = b\) krever et antall operasjoner som vokser omtrent som \((2/3) n^3\) med matrisestørrelsen.
For et triangulært system er beregningen mye enklere og vokser bare som \(n^2\) .
Forskjellen blir stor selv for moderate dimensjoner.

For å utnytte dette i Julia må vi markere matrisen som triangulær, for eksempel med UpperTriangular(U) eller LowerTriangular(L). Da velger Julia en langt raskere algoritme enn den som brukes for en vilkårlig matrise.

Eksemeplet nedenfor tester dette for en større matrise.

using LinearAlgebra, Random

n = 2000
Random.seed!(123)

A = randn(n, n)
b = randn(n)
U = triu(A)

println("\n--- Første kall (inkluderer kompilering) ---")
@time _ = A \ b;
@time _ = UpperTriangular(U) \ b;

println("\n--- Etter kompilering ---")
t1 = @elapsed for _ in 1:5
    _ = A \ b;
end

t2 = @elapsed for _ in 1:5
    _ = UpperTriangular(U) \ b;
end

println("Full matrise: $(round(t1/5*1e3, digits=2)) ms per løsing")
println("Triangulær matrise: $(round(t2/5*1e3, digits=2)) ms per løsing");

Oppgave: Gjør selv numeriske eksperiment og se hva slags skalering du får med n. 

Python: I python bruker man funksjonen solve_triangular fra scipy for å ta gevinst av denne farten

from scipy.linalg import solve_triangular
x = solve_triangular(U, b)

7.5.3 Systemer på formen \(A B x = c\)

Noen ganger møter vi systemer med produkt av matriser, for eksempel

\[ A B x = c. \]

Dersom både \(A\) og \(B\) er inverterbare, kan vi løse dette trinnvis:

\[ B x = y, \qquad A y = c. \]

Da får vi

\[ x = B^{-1} A^{-1} c. \]

Dette er nyttig i beregninger der én av faktorene har en kjent, enkel struktur
(for eksempel diagonal, ortogonal eller triangulær).

7.5.4 Ortogonal matrise

Om \(A\) er en ortogonal matrise, altså \(Q^\top Q = I\), får vi:

\[ Qx = b. \]

Løsningen er gitt direkte av

\[ x = Q^\top b. \]

Dette er numerisk stabilt og har en tydelig geometrisk tolkning:
\(Q\) endrer bare retning, ikke lengde — for eksempel en rotasjon eller refleksjon.
Å multiplisere med \(Q^\top\) roterer tilbake igjen.

Om \(A\) er en ortogonal matrise, altså \(Q^\top Q = I\), får vi

\[ Qx = b. \]

Løsningen er gitt direkte av

\[ x = Q^\top b. \]

Dette er numerisk stabilt og har en tydelig geometrisk tolkning:
\(Q\) endrer bare retning, ikke lengde — for eksempel en rotasjon eller refleksjon.
Å multiplisere med \(Q^\top\) roterer tilbake igjen.

Om vi bruker det faktum av kolonnene i \(Q\) er ortonormale,

\[ Q = [\,q_1 \; q_2 \; \dots \; q_n\,], \]

sår vi en interessant geometrisk tolking. Vi kan forstå løsningen som en projeksjon på hver av vektorene \(q_n\) i det ortonormale koordinatsystemet

\[ Qx =b = \sum_i q_i x_i \]

slik at \(x_i = q_i^\top b\) er vekten langs denne retningen i koordinatsystemet spent ut av kolonnene av \(Q\).

7.6 Ekstremalpunkter av flervariabel polynomfunksjoner

7.6.1 Gradient av polymomfunksjoner

En viktig klasse av funksjoner er polynomfunksjoner av flere variable, altså funksjoner som kan skrives som kombinasjoner av lineære og kvadratiske ledd i \(x\). Dette så vi blant annet at oppsto naturlig når vi lagde an Taylor-ekspansjon av funksjoner av flere variable.

Ofte ønsker man å derivere (det vil si å ta gradienten av) slike funksjoner, f.eks. for å finne maksiumum og minimum.

For er en lineær funksjon \[ f(x) = c^\top x, \] hvor \(c\) er en konstant vektor.

Gradienten blir da \[ \nabla f(x) = c. \]

Et naturlig neste steg er et kvadratisk former \[ g(x) = x^\top B x, \] der \(B\) er en symmetrisk matrise, da får vi \[ \nabla g(x) = 2 B x. \]

Bevis: Dette kan vi se ved å skrive ut uttrykket på elementform:

$$
g(x) = \sum_i \sum_j x_i b_{ij} x_j.
$$

Deriverer vi med hensyn på én komponent \(x_k\) får vi

$$
\frac{\partial g}{\partial x_k}
  = \sum_j b_{kj} x_j + \sum_i x_i b_{ik}.
$$

Siden \(B\) er symmetrisk (\(b_{ij}=b_{ji}\)), kan vi slå sammen leddene (og bytte on navnene på summeringsindeks):

\[ \frac{\partial g}{\partial x_k} = 2 \sum_j b_{kj} x_j. \]

Dermed får vi på matriseform

\[ \nabla g(x) = 2 B x. \]

Ta et 2×2-tilfelle med

\[ B = \begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix}, \qquad x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}. \]

Da er \[ g(x)=b_{11}x_1^2 + b_{12}x_1x_2 + b_{21}x_2x_1 + b_{22}x_2^2. \]

Ofte kombineres lineære og kvadratiske ledd:

\[ h(x) = x^\top B x + c^\top x + d, \]

som gir gradient

\[ \nabla h(x) = 2B x + c. \]

Høyere ordens polymer kan uttrykkes ofte som tensorer \(T_{ij...k}\) som er høyere ordens dimensjons-varianter av matriser. Ofte “flates” slike ut til matriser og vektorer, slik at man kan bruke linear algebra for å regne de effektivt.

7.6.2 Eksempel: Finne ekstremalpunkt for en to-variabel funksjon

La oss se på funksjonen

\[ f(x,y) = 3x^2 + 2xy + 2y^2 - 4x + y + 5. \]

Denne kan skrives på matriseform som

\[ f(u) = u^\top B\,u + c^\top u + d, \]

der

\[ u = \begin{bmatrix} x \\ y \end{bmatrix}, \qquad B = \begin{bmatrix} 3 & 1 \\ 1 & 2 \end{bmatrix}, \qquad c = \begin{bmatrix} -4 \\[2pt] 1 \end{bmatrix}, \qquad d = 5. \]

Vi ønsker nå finne ekstremala (maksimum eller minimum) av funksjonen.

Vi setter gradienten lik null:

\[ \nabla f(u) = 2Bu + c = 0. \]

som svarer til

\[ u= -\tfrac{1}{2} B^{-1} c. \]

Problemet er lite nok til å løses analytisk (med penn og papir), men vi løser det i Julia

# Først analytisk løsning
B = [3 1; 1 2]
c = [-4; 1]
d = 5

u = -0.5 * (B \ c)
println("Analytisk løsning:")
println(u)

# Evaluer funksjonen 
fval = u' * B * u + c' * u + d
println("Funksjonsverdi: ", fval)

Julia har også generiske numeriske løsere for generelle funksjoner. Dette minner om scipy.optimize (men drar nytte av auto-diff).

# numerisk løsning
using Optim

f(x) = x' * B * x + c' * x + d
result = optimize(f, [0.0, 0.0])
println("\nNumerisk løsning:")
println(result.minimizer)
println("Funksjonsverdi: ", result.minimum)

Normalt vet vi fra problemet om vi har et maksimum eller minimum. Men kanskje lyst å kontrollere eksplisitt funksjonen var faktisk er et minimum (etv. et maksimum) og visualisere formen på funksjonen. Dette kan vi gjøre med et contour-plot i 2d eller så kan vi plotte \(f(u - u_\mathrm{opt})\) som funksjon av \(\| u- u_\mathrm{opt}\|\) som fungerer for høyere dimensjoner, som vist nedenfor:

using LinearAlgebra, Plots

# --- Problem ---
B = [3 1; 1 2]
c = [-4; 1]
d = 5
f(x) = x' * B * x + c' * x + d

# Numerisk minimum
u_opt =  -0.5 * (B \ c)

# --- 1) Konturplott i (x,y) ---
x = range(-0.5, 2.0, length=200)
y = range(-2.0, 1.0, length=200)
z = [f([i, j]) for j in y, i in x]

contour(x, y, z;
    xlabel="x", ylabel="y", fill=true,
    title="f(x,y): konturer og minimum",
    aspect_ratio=1, cbar=false, label="")
scatter!([u_opt[1]], [u_opt[2]]; label="minimum")

# --- 2) Full grid → scatter: f vs ‖u - u_opt‖ ---
dists = []
vals  = []
for xi in x, yj in y
    u = [xi, yj]
    push!(dists, norm(u - u_opt))
    push!(vals,  f(u))
end

scatter(dists, vals;
    xlabel="‖u − u_opt‖",
    ylabel="f(u)",
    title="Fullt rutenett: f som funksjon av avstand",
    label="", ms=2)

contour contour

7.6.3 Optimering med bibetingelser

Vi ser nå på det samme problemet som tidligere, men legger til en bibetingelse (constraint). I mange anvendelser så optimerer vi funksjoner under slike vilkår, det kan f.eks. være

  • Størrelsen til en gjenstand. (Lengden på en stang)
  • Budsjetter som skal gå på i null
  • Masser og energi skal være bevart

Vi ønsker some tidligere å finne minimum av

\[ f(x,y) = 3x^2 + 2xy + 2y^2 - 4x + y + 5, \]
men nå under bibetingelsen

\[ x + y = 1. \]

eller under betingelsen

\[g(x,y) = x+y -1 =0\]

For akkurat dette problemet, så kan vi selvsagt bytte ut \(y\) med \(1-x\), og løse det som et en-dimensjonalt problem. Men for mer komplekse problem, som kan oppstå i praksis, er det ikke gitt at dette vil fungere. Om vi stedet optimerer funksjonen

\[ \mathcal{L}(x, y,\lambda) = f(x,y) + \lambda g(x,y) , \]

Så vil kravet \(g(x,y) = 0\) oppstå naturlig, når vi tar partiell-derivasjonen med hensyn på \(\lambda\). Dette systemet kan være enklere å løse. En slik måte å løse problemet på kalles å bruke en Langrange-multiplikator.

Som tidligere kan funksjonen kan som før skrives på matriseform:

\[ f(u) = u^\top B\,u + c^\top u + d, \]

der

\[ u = \begin{bmatrix} x \\ y \end{bmatrix}, \qquad B = \begin{bmatrix} 3 & 1 \\ 1 & 2 \end{bmatrix}, \qquad c = \begin{bmatrix} -4 \\[2pt] 1 \end{bmatrix}, \qquad d = 5. \]

Mens bibetingelsen skriver vi som \[ a^\top u - e =0 \] hvor \[ a = \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \qquad e = 1. \]

Da kan optimere funksjonen

\[ \mathcal{L}(u,\lambda) = u^\top B\,u + c^\top u + \lambda(a^\top u - e), \]

Vi setter så de deriverte lik null:

\[ \nabla_u \mathcal{L} = 2B\,u + c + a\,\lambda = 0, \qquad a^\top u = e. \]

Dette kan samles i ett lineært system:

\[ \begin{bmatrix} 2B & a \\ a^\top & 0 \end{bmatrix} \begin{bmatrix} u \\[2pt] \lambda \end{bmatrix} = \begin{bmatrix} -c \\[2pt] e \end{bmatrix}. \]


# ------------------------------------------------------------
# Optimalisering med en lineær bibetingelse
# f(x,y) = 3x^2 + 2xy + 2y^2 - 4x + y + 5
# Bibetingelse: x + y = 1
# ------------------------------------------------------------

# Problemdata
B = [3.0 1.0; 1.0 2.0]
c = [-4.0; 1.0]
d = 5.0

# Bibetingelse a' * u = d
a = [1.0; 1.0]
e = 1.0

# Sett opp systemet som skal løses
K   = [2B  a; a'  0.0]
rhs = [-c; e]

sol = K \ rhs
u   = sol[1:2]
lambda = sol[end]

# Funksjonsverdi i løsningen
f(u) = u' * B * u + c' * u + d

println("Løsning med bibetingelse: u* = ", u)
println("Lagrange-multiplikator:   lambda = ", lambda)
println("Funksjonsverdi f(u*) = ", f(u))
println("Sjekk bibetingelse a'u = ", a' * u, "  (skal være ", d, ")")

Dette kan vi plotte med

# --- Konturplot ---
x = range(-0.5, 2.0, length=200)
y = range(-2.0, 1.0, length=200)
fxy(x, y) = 3x^2 + 2x*y + 2y^2 - 4x + y + 5
z = [fxy(i, j) for j in y, i in x]

# Konturer
contour(x, y, z;
    fill=true, cbar=false, aspect_ratio=1,
    xlabel="x", ylabel="y",
    title="f(x,y) med bibetingelse x + y = 1",
    levels=20, label="")

# Linjen x + y = 1 (innenfor plottområdet)
x_line = range(-0.5, 2.0, length=200)
y_line = 1 .- x_line
mask = (y_line .>= minimum(y)) .& (y_line .<= maximum(y))
plot!(x_line[mask], y_line[mask], lw=2, label="x + y = 1")

# Marker løsningen
scatter!([u[1]], [u[2]]; color=:red, ms=6, label="løsning u*")

# Lagre eller vis
savefig("lagrange_contour.png")
display(current())

contour