computational methods

Question 1

what does it mean when dot = 0?

Answer 1

That vectors are orthogonal to each other.

Question 2

What is the goal of creating a polynomial?

Answer 2

To estimate the relationship between x and y by solving for unknowns so that a line (y) fits the data.

Question 3

What is the difference between a first and second degree polynomial?

Answer 3

1st degree ansatz is linear y = a0 + a1 * x

Second degree ansatz is y = a0 + a1 * x + a2 * x^2.

A second degree polynomial will give 3 columns in the matrix (coefficient, x, x^2) and 3 unknowns (a0, a1, a2).

Question 4

What does it mean that an equationsystem is overdetermined?

Answer 4

That you have more equations than unknowns to solve for. This usually happens when trying to fit a polynomial to a dataset. The data does not fit a straight line and there is no unique solution for the equation system. We will instead find the line that give the best estimate of the relationship of the data.

This is done by solving ATAv = ATy

Question 5

How do you create an equation system and fit a polynomial to a dataset?

Answer 5

Put the coefficients and x-values in one matrix.
- This matrix represents the equation system
a0 + a1 * x1 = y1
a0 + a1 * x2 = y2
ect.

If 2nd degree polynomial:
a0 + (a1 * x1) + (a2 * x2)

Put the unknowns in one vector
Put the y values in one vector.

This usually gives more equations than unknowns and we have to make the best estimate by solving ATAv = ATy.

Then use gaussian elimination (np.linalg.solve(ATA, ATy) to get the values for the unknowns and put them into the ansatz to get the polynomial.

Plug x-values into the polynomial. y = polynomial(x).

Question 6

What is the least squares method?
How do we solve for the residuals?

Answer 6

A way of finding the best estimate of a relationship between x and y where you minimize the residuals between the polynomial and the datapoints.

In other words to you find the relationship that is the best approximation of the real relationship.

r = y-vector - A-matrix * the values of the unknowns.

Question 7

What are norms?

Answer 7

The absolute size of vectors and matrices.

Question 8

What is a minimum norm, 2-norm and maximum norm?

Answer 8

minimum norm is 1-norm v1 + v2 +v3 …ect.
It does not take the “most efficient way”. Highest value of the columns in a matrix.

2-norm is called the euclidean distance and can be calculated with pythagoras. More efficient than minimum norm.

2-norm of v:
v^2 = x^2 + y^2.
vTv.

Maximum norm is the biggest value in the whole matrix (out of all columns and rows).

When calculating norms we always use absolute values.

Question 9

What is a condition number?

Answer 9

If you have an input error in your input data then it is going to be magnified in the output.

the condition number tells you how much it is going to multiply and is a sort of “warning sign” for how big your error will get.

We do not know the exact size of the difference between the real data and the input you used but we can calculate the norm of the difference.

Question 10

How do you calculate the relative error of your solution?

Answer 10

relative error = relative error in input data * condition number.

If your input is exact then input machine epsilon instead because in that case you have rounding error dominating.

Question 11

How can we make the condition number smaller?

Answer 11

By centering and scaling the data.

centering = subtract mean(x)
scaling = k( x - mean(x)) / sd(x)

ansats would be:
y = a0 + a1 * ( x - mean(x)) / sd(x)

Question 12

How do we get the condition number?

Answer 12

in python use np.linalg.cond(ATA).

Question 13

Why do we get rounding errors when doing computations in python?

Answer 13

A real number axis is infinite
computer number axis is discrete and finite and is representing infinity.

When you give python a number to store it is stored as the closest existing number on the computer number axis. This gives rounding errors every time we store a number.

Question 14

What are bits?

Answer 14

Numbers are stored in finite number of bits in binary form base 2 in computer memory.

Question 15

What is double precision?

Answer 15

64 bits.

sign (+/-) 1 bit.
exponent = 11bits
mantissa = 52 bits

Question 16

Define the terms precision and mantissa?

Answer 16

precision(p) = number of members in mantissa.

For example:
1.001 * 10^2

mantissa is 1.001, precision is 4 and base is 10.

Question 17

If we have mantissa 1.00 and 1.25 with base = 2, upper limit = 2 and lower limit = 0 of the exponent, what numbers can be stored? What is the precision?

Answer 17

1.00 * 2^0, 1.00 * 2^1, 1.00 * 2^2

1.25 * 2^0, 1.25 * 2^1, 1.25*2^2

precision = 3

Question 18

What is the difference between binary numbers and decimal numbers?

Answer 18

binary has base 2
decimal has base 10

Question 19

How can you convert a binary number to a decimal number?

Answer 19

Ex binary number 1.001

decimal = 12^0 + 02^-1 + 02^-2 + 12^-3

decimal = 1.125

Take the number in the mantissa multiplied by base with exponent like the index position of number.

Question 20

How do you calculate absolute error and relative error? Do they stay the same across the entire computer number axis?

Answer 20

abs = a - b where a is the true value

rel = a - b / a where a is the true value

abs error get bigger with higher numbers but relative error is the same across the entire axis = machine epsilon.

Question 21

What is machine epsilon?

Answer 21

The magnitude of the relative round off error and it is usually 10^-15.

Question 22

How can we check if two values can be considered equal?

Answer 22

Look to see if the relative difference between the stored value and the truth is smaller than machine epsilon.

Question 23

What is machine epsilon?

Answer 23

The magnitude of the relative roundoff error. 10^-16

Question 24

What are ODEs?

Answer 24

Models used to describe change, often in time. The model is written as derivatives.

The solution y(t) is unknown and cannot be found. But we can use numerical methods like Euler’s and Heun’s to simulate the changes of y with time and get an approximation of y(t) by using y’(t) which is the rate of change.

Question 25

When using an ODE to get an approximation of y(t), what is h?

Answer 25

h is the distance between the time steps.
Using smaller h will improve accuracy but will give more computations since we have to calculate the solution at every timepoint.

Because of this we want to use numerical methods that do more calculations per step which gives higher accuracy without having to lower h.

Question 26

How do we use ODEs to solve for y(t) in python?

Answer 26

The derivative is given as a function which is then used in a solver as input. If you have more than one derivative and solution then give the derivatives as a vector.

Def (t, y, constants)
X, S = y
constants = a, b, c
yt = [x’(t), s’(t)]

return yt

solver is solve_ivp.

Question 27

What will be the output of solve_ivp?

Answer 27

A time vector with the timepoints.

A matrix with the solutions from the different derivatives as rows and time points where the solutions were stored as columns.

If you have 9 derivatives you will get 9 rows in the matrix.

Question 28

What is the difference between Euler’s and Heun’s method?

Answer 28

To get the y(t) in Euler’s method we calculate the tangent in only f(yi, ti).
- gives one calculation per timestep.

To get y(t) in Heun’s we find the tanget in both f(yi, ti) and f(yi+1 , ti+1) and take the average of them.
- gives double the amount of calculations per time step but giver higher accuracy and we can use bigger h to keep error tolerance.

Question 29

What is discretization and discretization error?

Answer 29

Since we can’t compute infinity our approximations need to represent infinity. This gives many approximations.

Discretization is approximations in discrete points.

Discretization error is the error that occur because of the discretization and is dependent on h. If we use a smaller h the error decreases.

Question 30

Discretization error can be divided into two parts, what are they?

Answer 30

Local and global error.

Local error is the error introduced in one step.

Global error is the total error in every step.

Question 31

What is Taylor’s expansion?

Answer 31

Used to give a mathematical explanation for the discretization error.

Question 32

What does Taylor expansion say about the discretization error for Euler’s vs Heun’s vs classical range kutta?

Answer 32

Eulers:
error in one step = c1 *h^2 + c2 * h^3….ect.

c1*h^2 dominates since it is the biggest and then the error gets smaller and smaller.

the global error becomes n(c*h) where n is number of time steps (tend - t0 / h).

Heun’s:
total error is c * h^2

Classical range kutta:
total error is c*h^4

Question 33

What is order of accuracy?

Answer 33

The discretization error decreases at different rates depending on method. How much the error decreases with every step is called the order of accuracy where higher order indicates a better method because you can have bigger h and still keep a certain error tolerance.
There are more computations per step but it is still going to be faster than many time steps.

Eulers: first order, needs small h
Heun’s: second order, can have bigger h
classical range kutta: fourth order, can have even bigger h.

Question 34

What is the big O notation?

Answer 34

A way to write the order o a method.

Eulers O(h)
Heun’s O(h^2)

Question 35

What is adaptivity?

Answer 35

We need adaptivity to keep accuracy and it means that the value of h changes with how big the derivatives are because the discretization error gets bigger with bigger derivatives.

If the changes are big then h gets smaller to stay within the error tolerance.

Question 36

How does the adaptivity algorithms usually work?

Answer 36

It calculates a new value and the local error.

If the local error is bigger than the given tolerance then reduce h and do the calculation again.

Question 37

How do you calculate the new error when you reduce h by a factor of 10 for

Euler’s
Heun’s
Cl.RK

Explain why

Answer 37

When you want to calculate the new error you divide the original error by the factor you are reducing h by to the power of the order of the method.

This because the discretization error is dependent on h and is proportional to the order of the methods.

e1/e2 = (h1/h2)^p where p is the order of the method.

Question 38

Can you solve higher derivatives with solve_ivp?

Answer 38

Yes but you have to rewrite the ODEs so that you don’t have a derivative in the derivative you give to solve_ivp.

You rewrite the derivative so that you have:
y’(t) = f(t,y)
y’‘(t) = f(t, u1, u2)

then you can use solve_ivp

Question 39

What is the difference between Monte Carlo and numerical methods to solve ODEs?

Answer 39

Monte Carlo methods are stochastic and ODEs are deterministic.

Monte Carlo methods will give different results to the same thing each time because it is based on randomness. ODEs will give the same result each time because the solution is predictable and determined given a certain input data.

Question 40

When is it better to use monte carlo methods vs deterministic methods?

Answer 40

If you want to look at the change of something over time and you have very few datapoints then use a monte carlo simulation. This because the change may happen at t(1) or t(8) (almost random) and the deterministic method will most likely be off. When you look microscopically at the individual molecules.

If you have many datapoints then use a deterministic method because it is more accurate for simulations of many datapoints. Consider however that with many trials, monte carlo methods give results close to the deterministic but it is not efficient for too many datapoints.

You may also consider that deterministic methods gives decimal results which in reality might be impossible.

Question 41

Explain and exemplify deterministic methods and models.

Answer 41

A deterministic method (such as Euler’s Heun’s) is a technique or algorithm that, given the same initial conditions and inputs, will produce the same result every time it is executed. There is no randomness involved in the process. The time steps and equations are fixed inputs.

A deterministic model is a mathematical representation of a system or a process where the output is entirely determined by the initial conditions and input parameters. For example ODEs

Question 42

Why does monte carlo methods get higher accuracy with multiple tries?

Answer 42

Becuse of the central limit theorem. The standard deviation will decrease with higher number of tries. The slope here is the same as the order of accuracy.

This is the basis of MC-methods.

Question 43

What can monte carlo methods be used for?

Answer 43

They do the same thing as ODEs, obtain numerical results that approximate results at certain times but it uses probability and randomness to do so.
It calculates the probability of reaction/change and going from one state to another each trial.

This is called a stochastic method.

Question 44

What is a propensity function in a monte carlo method?

Answer 44

The function for the probability of moving from one state to another.

Question 45

What is the general process of monte carlo methods?

Answer 45

Input (N), perform multiple MC simulations and get the results, take the average of the results.

The more trials you do, the higher the accuracy will be.The bell shape of the results gets narrower and gathers around the truth.

Question 46

How does the standard deviation decrease with increasing number of trials in MC simulations?

Answer 46

It decreases according to the central limit theorem.

1 / sqrt(N).

The slope of this curve can be compared to the order of accuracy for methods used to solve ODEs.

the convergence of monte carlo methods are slower than for ex. Euler’s.

Question 47

Explain what a Markov process is.

Answer 47

It is a stochastic process meaning that it uses randomness. It also satisfies the Markov property that says “ One can make predictions for the future of the process based solely on its present state”

Question 48

What is the SSA algorithm?

Answer 48

Stochastic simulation algorithm/Gillespie algorithm.

It is a Monte Carlo method and a discrete stochastic Markov chain simulation used to simulate change between states using probabilities and randomness.

The time steps and the reactions are chosen based on randomness and probabilities.

Question 49

What input data does the SSA algorithm need?

Answer 49

propencity vector - has the propencities that become probabilities when divided by a0. Values will change with every iteration when states changes.

state vector - the states to calculate the probabilities for moving through. Changes every iteration.

Stoichiometry matrix - How something changes in a reaction. Does it increase or decrease in the reaction?

Question 50

When we use two different numerical methods to solve an ODE with different order of accuracy, why don’t we get two different results since they differ in accuracy?

Answer 50

Because most methods are adaptive and they will change h to keep an error tolerance. The one with lower order will be much slower however because it will need very small time steps and do the calculations many more times.

Question 51

What are the different steps of the SSA algorithm?

Answer 51

1. Find time to next reaction.
2. Find next reaction. Using the CDF.
3. Update state vector.
4. Update the time.

Question 52

If the starting values are really high, why would it not be good to use a monte carlo method?

Answer 52

Because really high starting values will give really high a0, which in turn gives very small timesteps. It would have to take many steps before anything happens and it can get inefficient for a really large number of units.

Question 53

What is a0 in the SSA algorithm and what happens if it is large vs small?

Answer 53

a0 gives the steepness of the exponential distribution curve that we sample random time steps from.

If it is large then the timesteps will be really small and if a0 is small then the timestep will be bigger.

Question 54

Why is it bad to have to small a0?

Answer 54

Because smaller a0 will give large time steps which in turn mean that it will take longer until any reaction. You could also step outside of the time range and nothing happens at all.

Question 55

What is a PDF and CDF?

Answer 55

The PDF is the probability distribution function and describes the probabilities for for example different reactions to happen.

The CDF is the cumulative distribution function and gives the probability that a random variable is less than or equal to a certain value.

It is used to go from a uniform distribution to any other using the inverse transform sampling. This is what we do in the SSA algorithm to get the next reaction.

Question 56

How does the inverse transform sampling differ in discrete and continuous cases?

Answer 56

Is discrete cases we sum the probabilities and in continuous cases we integrate them.

Question 57

Describe how you get the probability that X is less or equal to x in a discrete case.

Answer 57

Sum the probabilities and then use the inverse transform sampling.

Choose a number from the uniform distribution on [0,1] and plug it in to the y-axis of the CDF and get a number from another distribution.

The CDF will differ depending on what distribution you want.

Question 58

Describe how you get the probability that X is less or equal to x in a continuous case.

Answer 58

The same way as in a discrete case but the sum become integration because the “jumps” are infinite.

The inverse transform sampling is:

u = 1 - e^-a0T
Will give:
T = - 1/a0 * ln (1-u)

Plug in u and get T where u is the random number from the uniform distribution.

Question 59

Give examples of deterministic and stochastic models and methods that solves them?

Answer 59

Deterministic model: ODE
Stochastic model: Brownian motion

Deterministic method: Heun’s, Euler’s
Stochastic method: SSA, monte carlo

Question 60

What is the difference between a model and a method?

Answer 60

A model is the mathematical description of how something works. Like ODEs or brownian motion while the methods are how we solve the problems like Heun’s or SSA.

Question 61

You are running a monte carlo simulation and get a certain confidence interval. How can you reduce it by factor 1/4?

Answer 61

Increase number of trials by 4^2 where.

This is because the sd decreases by 1 / sqrt(N) because of the central limit theorem.

To reduce it by 4 you take 1/4 * 1 / sqrt(N) which gives factor 4^2.

Question 62

What is min∣∣b−Az∣∣
2?
​

Answer 62

Solving for min∣∣b−Az∣∣
2 means that your solving for the vector z that minimizes the euclidean norm of the residual vector b-Az.

It is the least squares method.

Question 63

How would you calculate the new h if your error is reduced by a factor 10 for

Euler’s
Heun’s
Cl.RK

Answer 63

If you want to calculate how h changes with reduced error you instead divide the original h by the factor the error is reduced by squared.

Question 64

Explain and exemplify stochastic models and methods.

Answer 64

A stochastic method (such as SSA, monte carlo methods) use randomness and probability to get the results. Due to the randomness multiple trials will not give the same results. Both time to next reaction and reaction/equations are chosen randomly in these methods.

A stochastic model (such as brownian motion) is mathematical representation of a system or process that incorporates randomness or uncertainty.

Question 65

Explain brownian motion

Answer 65

Brownian motion is a stochastic continuous Markov process that describes the random movements of particles due to them crashing together.

The applications are however widespread and the model is commonly used in monte carlo methods.

In brownian motion each time step is normally distributed.

Question 66

Imagine we have a curve of the discretization error for a method. What in the curve describes:

• Discretization error
• Rounding error
• Order of accuracy
• Machine epsilon

Answer 66

Discretization error is represents by the curve until it hits machine epsilon 10^-15, because it cannot get smaller than that. The curve after (the positive turn) describes how the rounding errors start to dominate.

The slope is order of accuracy

Question 67

Is SSA a monte carlo method?

Answer 67

It is treated as such but actually the SSA algorithm does not do the code many times so strictly it is not a monte carlo method.

Question 68

If you know the time it takes to calculate one time step for Euler’s method, how can you find out how long one step takes for Heun’s and Cl.RK?

Answer 68

Heun’s has double amount of calculations per step and takes twice as long as Euler’s to calculate one step. However the total time will be shorter because the number of timesteps are fewer for Heun’s than Euler’s.

Cl.RK has 4 times as many calculations as Euler’s and each timestep takes for times as long as Euler’s but the total time will be much shorter due to fewer steps.

Question 69

How do we calculate the number of time steps in a numerical method?

Answer 69

time interval / h

Question 70

What is ||truth - reality || divided by ||truth||?

Answer 70

The relative error of a solution.

Question 71

Explain what Range-Kutta method of order 5(4) is.

What different errors do we have in these methods, which one is controlled?

Answer 71

A numerical method Range-Kutta has two computations running in parallel with each other.

One where the discretization error is dependent on h^4 and the other dependent on h^5.

Accuracy is based on order 4.

In computations like these we have discretization error and rounding errors, we are controlling the discretization error.

Question 72

Does “two methods are running simultaneously” mean that the computational cost is doubled (compared with just
running one “pure” method?

Answer 72

The computational cost is not doubled.

The methods are constructed so that the ki’s are used in both methods and does not need to be calculate for both methods. There is only a small overhead this way.

Question 73

What input does the SSA algorithm need?

Answer 73

propencity vector
stochiometry vector
y0 (state vector)
time interval
parameters