Programmering

SNDtrapes.py

# -*- coding: utf-8 -*-
"""
Created on Sun Nov 27 10:34:39 2022
 
@author: Sebas
"""
 
import numpy as np
 
def SND(x): #Dette er funksjonen vår, standard normalfordeling
    return 1/np.sqrt(2*np.pi) * np.exp(-(x**2)/2)
 
def trapesmetoden(f, a, b, N=1000): 
    #Dette er en funksjon for å ta trapesmetoden mellom a og b for en hvilken som helst funksjon "f" som gis som argument
    h = (b-a)/N #steglengde mellom a og b med N steg
    x = np.linspace(a,b,N+1) #punktene til x mellom a og b
 
    res = 0 #tellevariabel
    for i in range(N):
        res += (f(x[i])+f(x[i+1]))/2 * h #trapesmetodens definisjon
 
    return res
 
print(trapesmetoden(SND, 0, 1)) 
#Print standard normalfordelings sannsynlighet for et resultat mellom 0 og 1 (innenfor ett standardavvik over gjennomsnittet)
#Forventet resultat er ~34.1% (altså: 0.341...)

ellipseOmkrets.py

# -*- coding: utf-8 -*-
"""
Created on Tue Nov 22 08:48:23 2022
 
@author: Sebas
"""
 
import numpy as np
import matplotlib.pyplot as plt
 
a = 2 # halvakse 1
b = 1 # halvakse 2
N = 1000 #antall steg i trapesmetoden
Dx = a/N #steglengde
 
def df(x): #Den deriverte av en parametrisering y=f(x) for ellipsen.
    return -b/(a**2) * (1/(np.sqrt(1-(x**2)/(a**2)))) * x
 
def phi(x): #Buelengde-funksjonen som skal integreres
    return np.sqrt(1+df(x)**2)
 
def trapesmetoden(): #trapesmetoden
    tellevariabelen_min = 0
    for i in range(N-1):
        x_i = 0 + Dx*i #0 er startpunktet vi integrerer fra
        tellevariabelen_min += (phi(x_i + Dx) + phi(x_i))*Dx/2 #Trapesmetode-definisjonen
    return tellevariabelen_min
 
#Vi fant bare buelengden av en fjerdedel av ellipsen, så vi ganger med 4
print(trapesmetoden()*4)

NumeriskDifLikning.py

# -*- coding: utf-8 -*-
"""
Created on Tue Nov  8 11:25:09 2022
 
@author: Ylva
"""
 
import numpy as np
 
a = 0.00001 #start-x-verdi
b = np.pi #slutt-x-verdi
N = 10000 #antall x-punkter
 
dX = (b - a)/N #lengden mellom punktene
 
def sinc(x):
    return np.sin(x)/(x) #sinc kan også defineres som sin(pi*x)/(pi*x). Livet er trist
 
def trapes(f): #MERK: denne funksjonen tar inn en funksjon! KULT!
    area = 0
    for i in range(N):
        area += 1/2 * (f(a + dX * i) + f(a + dX * (i + 1))) * dX
        #Her går det ann å skrive x = a + dX*i og legge det inn 
        #i funksjonen for å være enda tydeligere.
    return area
 
print(trapes(sinc))
 

NumeriskDifLikning.py

# -*- coding: utf-8 -*-
"""
Spyder Editor
 
This is a temporary script file.
"""
 
import numpy as np
 
def f(x):
    return np.exp(- x**2) #funksjonen vår
 
a = 0 #start-x-verdi
b = 1 #slutt-x-verdi
n = 1000 #antall punkter på x
 
res = 0 #variabel som lagrer summen "res" for "resultat"
 
dx = (b-a)/n #Delta x
x = np.zeros(n) #Våre x-verdier
#alternativ: x = np.linspace(a,b,n) #Trenger da ikke x i for-løkke
 
 
for i in range(0,n):
    x[i] = a + i*dx #MERK: starter i a
    res += f(x[i]) #Legger til rektangelets høyde
    #Skal egentlig også gange med dx, men tar det senere
    #for å spare prosessorkraft
 
res *= dx #ganger med dx
print("Riemann:",res)
 
 
#Trapesmetoden på Seb-måten:
res = 0
dx = (b-a)/n #Delta x
x = np.linspace(a,b,n)
for i in range(0,n-1):
    #antar vi brukte linspace tidligere
    res += (f(x[i])+f(x[i+1]))/2
res *= dx
print("trapes:",res)

NumeriskDifLikning.py

import numpy as np
from tqdm import tqdm #Kun for å få progress bar, fordi det ser kult ut
import matplotlib.pyplot as plt
 
# g(x) er det som står på høyre side av likningen når x' står alene på venstre side.
def g(x):
    return -(1/x**2)*((np.e**(1/x))/((np.e**(1/x)-1)**2))*(x-100) #Siste tall (100) er T_k.
 
#Vi definerer den derriverte av g(x) siden den brukes i noen av metodene.
def gDeriv(x):
    return (g(x+0.0001)-g(x-0.0001))/0.0002
 
 
x0 = 1  # x(0)
maxT = 10  #hvor stor definisjonsmengden til grafen vi skal plotte er
steps = 10  #hvor mange skritt vi tar, høyere verdi = bedre resultat men saktere
 
h = maxT/steps #Lengden av hvert skritt
 
def eulersMethod():
    x = [x0] #Lager en liste med y-verdier til grafen til x(t), første element er initialverdien x0
    ts = [0]  #lager en liste med x-verdier
 
    for n in tqdm(range(steps)):
        ts.append((n+1)*h) #Listen med x-verdier starter på 0, og for hver iterasjon legger vi på (n+1)*h, 
                           # slik at den ender opp som [0, 1h, 2h, 3h, 4h ... ]
        x.append(x[n] + h*g(x[n])) #Her bruker vi eulers eksplisitte metode til å finne hva neste x er, og legger
                                   # den til i lista som heter x. 
 
    return ts, x #Her returneres de to listene slik at de kan plottes (se bunnen av koden)
 
 
def eulerImplicitMethod():
    x = [x0]
    ts = [0]
 
    for n in tqdm(range(steps)):
        ts.append(n*h+h)
 
        #Newtons metode
        #Her brukes newtons metode til å finne x_(n+1), ved å løse likningen x_(n+1) = x_n + h*g(x_(n+1)).
        #1: Skriv om likningen til 0 = -x_(n+1) + x_n + h*g(x_(n+1)), det som står på høyre side er nå f(x) funksjonen vi bruker i newtons metode.
        #   altså, f(x) = -x_(n+1) + x_n + h*g(x_(n+1))
        #2: Finn f'(x) ved å derrivere det over, da får vi f'(x) = h*g'(x_(n+1)) - 1
        #3: Hvis vi nå setter inn f(x) og f'(x) i formelen for newtons metode, får vi at det man må gjøre hver
        #   iterasjon er: x_(n+1) = x_(n+1) - (-x_(n+1) + x_n + h*g(x_(n+1))) / (h*g'(x_(n+1)) - 1)
 
        nextX = x[n] #Denne variabelen er x_(n+1) (det er ikke lov å bruke x_(n+1) som variabelnavn i python... )
        for i in range(100):
            lastX = nextX #Dette er kun for å sjekke om newtons metode har konvergert
 
            nextX = nextX - ((-nextX+x[n]+h*g(nextX))/(h*gDeriv(nextX)-1)) #Her gjør vi iterasjonen
 
            if abs(nextX - lastX) < 0.000000001:
                break  #Dersom den har konvergert avslutter vi loopen (for å ikke bruke unødvendig lang tid)
 
        x.append(nextX) #Nå når vi har funnet x_(n+1) med newtons metode legger vi den til i lista x
 
    return ts, x
 
 
def trapesMethod():
    x = [x0]
    ts = [0]
 
    for n in tqdm(range(steps)):
        ts.append(n*h+h)
 
        #Newtons Metode
        nextX = x[n]
        for i in range(100):
            lastX = nextX
            nextX = nextX - ((-nextX+x[n]+(h/2)*g(x[n])+(h/2)*g(nextX))/((h/2)*gDeriv(nextX)-1))
            if abs(nextX - lastX) < 0.000000001:
                break
 
        x.append(nextX)
 
    return ts, x
 
 
def heunsMethod():
    x = [x0]
    ts = [0]
 
    for n in tqdm(range(steps)):
        ts.append(n*h+h)
        x.append(x[n] + (h/2)*(g(x[n]) + g(x[n] + h*g(x[n]))))
 
    return ts, x
 
 
#Her plottes det en løsning med eulers metode
xs, ys = heunsMethod()
plt.plot(xs, ys, color="blue")
 
 
#Her plottes det en løsning med trapesmetoden
xs, ys = trapesMethod()
plt.plot(xs, ys, color="red")
 
 
#Her endrer vi steps til et veldig høyt tall og plotter løsning med heuns metode for å få en "fasit" på hvordan 
#grafen skal se ut, slik at vi kan se hvilken av eulers metode og trapesmetoden som er best.
steps = 1000000
h = maxT/steps
xs, ys = heunsMethod()
plt.plot(xs, ys, color="green")
 
 
plt.show()

Oppgave 13:

from math import exp
import matplotlib.pyplot as plt
 
 
def g(T):
    return 2 - T
 
 
T_0 = 1 / 10
h = 0.01
n = 1000
 
T_eksplisitt = [T_0]
t = [0]
 
x = [0, 1, 2, 3]
x[2]
 
 
for i in range(n):
    # Euler eksplisitt
    T_eksplisitt.append(T_eksplisitt[i] + h * g(T_eksplisitt[i]))
    t.append((i + 1) * h)
 
plt.plot(t, T_eksplisitt)
plt.show()

Oppgave 14:

from math import exp
import matplotlib.pyplot as plt
 
 
def g(T):
    return (T_k - T) * 1 / T ** 2 * exp(1 / T) / (exp(1 / T) - 1) ** 2
 
 
T_k = 2
T_0 = 1 / 10
 
n = 1000
h = 0.01
 
T_eksplisitt = [T_0]
T_heuns = [T_0]
x = [0]
 
for i in range(n):
    # Euler eksplisitt
    T_eksplisitt.append(T_eksplisitt[i] + h * g(T_eksplisitt[i]))
    # Heuns metode
    T_heuns.append(T_heuns[i] + h * ((g(T_heuns[i]) + g(T_heuns[i] + h * g(T_heuns[i]))) / 2))
 
    x.append((i + 1) * h)
 
plt.plot(x, T_eksplisitt)
plt.plot(x, T_heuns)
plt.show()

Øl-graf: https://www.desmos.com/calculator/alhga7sjdb

    #%% Den tingen som står på toppen
# -*- coding: utf-8 -*-
"""
Created on Thu Sep  8 15:18:33 2022
 
@author: Sebastian Westerlund Kalland
"""
 
#%% Imports
import numpy as np
import matplotlib.pyplot as plt
 
#%% Felles variabler
N = 10
L = 1
h = L/N
a = 2
t = np.linspace(0,L,N+1)
 
#%% Vanlig graf
def x(t):
    return np.exp(a*t)
 
#%% dx
def dx(x):
    # x derivert er lik a ganger x. Altså, denne funksjonen representerer x-derivert ved en verdi x
    return a*x
 
#%% Euler eksplisitt
x_n_exp = np.zeros(N+1)
x_n_exp[0] = 1
 
for i in range(N):
    x_n_exp[i+1] = x_n_exp[i] + dx(x_n_exp[i])*h
 
#%% Euler implisitt
x_n_imp = np.zeros(N+1)
x_n_imp[0] = 1
 
for i in range(N):
    #x_n_imp[i+1] = x_n_imp[i] + h*dx(x_n_imp[i+1]) 
    #LØS FOR x_n_imp[n+1]
    #Her løst for dx = a*x:
    x_n_imp[i+1] = x_n_imp[i]/(1-a*h)
 
#%% Trapesmetoden
x_n_tra = np.zeros(N+1)
x_n_tra[0] = 1
 
for i in range(N):
    #x_n_tra[i+1] = x_n_tra[i] + h*(dx(x_n_tra[i]) + dx(x_n_tra[i+1]))/2 
    #LØS FOR x_n_tra[n+1]
    #Her løst for dx = a*x:
    x_n_tra[i+1] = x_n_tra[i]*(1+h*a/2)/(1-a*h/2)
 
#%% Midtpunktmetoden
x_n_midt = np.zeros(N+1)
x_n_midt[0] = 1
 
for i in range(N):
    #x_n_midt[i+1] = x_n_midt[i] + h*dx((x_n_midt[i] + x_n_midt[i+1])/2)    
    #LØS FOR x_n_midt[n+1]
    #Her løst for dx = a*x:
    x_n_midt[i+1]= x_n_midt[i]*(1+h*a/2)/(1-h*a/2)
    #I akkurat dette tilfellet blir midtpunktmetoden og trapesmetoden like. Kan du se hvorfor?
 
#%% Heuns metode
 
x_n_H = np.zeros(N+1)
x_n_H[0] = 1
 
for i in range(N):
    estimat = dx(x_n_H[i] + dx(x_n_H[i])*h)
    x_n_H[i+1] = x_n_H[i] + h*(dx(x_n_H[i]) + estimat)/2
 
#%% Masse figurting, bare ignorer
plt.figure(dpi=300)
plt.rc('axes',labelsize=12)
plt.xlabel('$t$')
plt.ylabel('$x$')
plt.tick_params(axis="both", length=3)
plt.minorticks_on()
plt.grid(which="minor", axis="both", linestyle=":")
plt.grid()
 
#%% Plotting
plt.plot(t, x(t), label='$x(t)=e^{at}$')
 
plt.plot(t, x_n_exp, label='Euler eksplisitt for $\dot{x}=ax$')
 
plt.plot(t, x_n_imp, label='Euler implisitt...')
 
plt.plot(t, x_n_tra, label='Trapesmetoden...')
 
plt.plot(t, x_n_midt, label='Midtpunktmetoden...')
 
plt.plot(t, x_n_H, label='Heuns metode...')
 
#%% liten ekstra plotte-ting
plt.legend()

    # -*- coding: utf-8 -*-
"""
Created on Wed Sep  7 17:41:55 2022
 
@author: Sebastian Westerlund Kalland og Truls
"""
 
import numpy as np 
def newtonsmetode(func, der_func, x=1, n=10, d=0):   
    '''
    Parameters
    ----------
    func : String
        0 = "func". This is a function expressed as string using x as its input. Example: "x - np.cos(x)"
    der_func: String
        Derivative of f.
    x : float
        Starting value. Default 1
    n : integer
        Maximum number of iterations. Default 10.
    d : float
        Accepted error margin before quitting the iteration prematurely. Default 0, meaning no errors allowed.
 
    Returns
    -------
    x :
        Value after n iterations of Newton's method
    '''
    f = lambda x: eval(func)
    df = lambda x: eval(der_func)
 
    for i in range(n):
        x_pre = x
        x = x - (f(x)/df(x))
        if abs(x_pre-x)<d:
            print("x =",x,"\tafter",i,"iteration, with an error margin of",d)
            return x
        if abs(x_pre-x)==0:
            print("x =",x,"\tafter",i,"iterations, the computer cannot represent x more accurately.")
            return x
 
    print("x =",x,"\tafter",n,"iterations.")
    return x
 
newtonsmetode("x-np.cos(x)", "np.sin(x)+1", 1, 20, 1e-6)
newtonsmetode("x-np.cos(x)", "np.sin(x)+1", 1, 3)
newtonsmetode("x-np.cos(x)", "np.sin(x)+1")

    # -*- coding: utf-8 -*-
"""
Created on Tue Sep  6 11:22:31 2022
 
@author: Sebastian Westerlund Kalland
"""
 
import numpy as np
 
x = 1
for i in range(100):
    forrige_x = x
    x = np.sqrt(np.cos(x))
    print(i, x)
    if abs(forrige_x-x) < 0.000001:
        break
 
def f(x):
    return np.cos(x) - x**2
 
def df(x):
    return -np.sin(x) - 2*x
 
x = 1
for i in range(10):
    forrige_x = x
    x = x - f(x)/df(x)
    print(i, x)  
    if abs(forrige_x-x) < 0.000001:
        break
 
#while-eksempel for de som er litt mer spenstig
x = 1
forrige_x = 2
while abs(forrige_x-x) >= 0.000001:
    forrige_x = x
    x = x - f(x)/df(x)
    print(x)

import numpy as np
from math import factorial, pi
from sympy import diff, Symbol, lambdify, sin, cos
from matplotlib import pyplot as plt
 
# Sets variable
x = Symbol('x')
 
 
def taylor(f, n: int, a: int, x_values: list[float]):
    """
    ´n´: Degree of taylor polynomial
    ´a´: X-value at which our approximation is the most exact
    """
    num_values = len(x_values)
    # Fills in values for f(a)
    y_values = np.ones(num_values) * lambdify(x, f, "numpy")(np.array([a]))[0]
    # Differentiates function n times and adds it all up
    for i in range(n):
        f = f.diff(x)
        f_prime_a = lambdify(x, f, "numpy")([a])[0]
        y_values += f_prime_a / factorial(i + 1) * lambdify(x, (x - a) ** (i + 1), "numpy")(x_values)
    return y_values
 
 
def main():
    # Defines function we want to approximate
    f = cos(x)
 
    x_values = np.linspace(-2 * pi, 2 * pi, 100)
    y_values = taylor(f, 4, 0, x_values)
    plt.plot(x_values, lambdify(x, f, "numpy")(x_values))
    plt.plot(x_values, y_values)
    plt.ylim([-5, 5])
    plt.show()
 
 
if __name__ == '__main__':
    main()