LAB EXERCISE 6

Aim: Runge Kutta’s 2nd order Program in R

Theory

The Runge-Kutta method finds an approximate value of y for a given x. Only first-order
ordinary differential equations can be solved by using the Runge-Kutta 2nd-order method.
Below is the formula used to compute the next value yn+1 from the previous value yn.
Therefore:

yn+1 = value of y at (x = n + 1)

yn = value of y at (x = n)

where

0 ? n ? (x - x0)/h

h is step height

xn+1 = x0 + h

The essential formula to compute the value of y(n+1):

K1 = h * f(x,y)

K2 = h * f(x/2, y/2) or K1/2

yn+1 = yn + K2 + (h3)

The formula basically computes the next value yn+1 using current yn plus the weighted
average of two increments:
· K1 is the increment based on the slope at the beginning of the interval, using
y.
· K2 is the increment based on the slope at the midpoint of the interval, using
(y + h*K1/2).

The method is a second-order method, meaning that the local truncation error is in the order
of O(h3), while the total accumulated error is order of O(h4).
Code

f=function(x,y){

y*sin(x)

RK2 = function(f, xn, yn, h, N, n=0){ #function to calculate the differential

df = NULL

M = f(xn,yn)

while(n < N){

df = rbind(df, [Link](n=n, Xn = xn, Yn = yn, M = M, hM = h*M))

#to create the data frame of the values

K1 = M

K2 = f(xn+h, yn+K1)

M = (K1 + K2)/2 #value of differential

xn = xn + h #updating the values

yn = yn + h*M

n=n+1

print(df)

df = RK2(f, -3, 1, 0.1, 5) #calling the function

Conclusion

Implemented Runge Kutta 2nd order rule in R language, along with its theory and
mathematical formulae.
 LAB EXERCISE 7

Aim: Runge Kutta’s 4th order Program in R

Theory

Given the following inputs,

Thus we are given below.

The task is to find the value of the unknown function y at a given point x.
The Runge-Kutta method finds the approximate value of y for a given x. Only first-order
ordinary differential equations can be solved by using the Runge Kutta 4th order method.
Below is the formula used to compute next value yn+1 from previous value yn. The value of n
are 0, 1, 2, 3, ….(x – x0)/h. Here h is step height and xn+1 = x0 + h
Lower step size means more accuracy.

The formula basically computes next value yn+1 using current yn plus weighted average of
four increments.
· k1 is the increment based on the slope at the beginning of the interval, using y
· k2 is the increment based on the slope at the midpoint of the interval, using y + hk1/2.
· k3 is again the increment based on the slope at the midpoint, using y + hk2/2.
· k4 is the increment based on the slope at the end of the interval, using y + hk3.

The method is a fourth-order method, meaning that the local truncation error is on the order
of O(h5), while the total accumulated error is order O(h4).

Code

f=function(x,y){

y*sin(x)

RK4 = function(f, xn, yn, h, N, n=0){ #function to calculate the differential

df = NULL

M = f(xn,yn)

while(n < N){

#to create the data frame of the values

K1 = M

K2 = f(xn+h/2, yn+K1/2)

K3 = f(xn+h/2, yn+K2/2)

K4 = f(xn+h, yn+K3)

M = (K1 + 2*K2 + 2*K3 + K4)/6 #value of differential

xn = xn + h #updating the values

yn = yn + h*M

n=n+1

print(df)

df = RK4(f, -3, 1, 0.1, 5)

Output

Conclusion

Implemented Runge Kutta 4th order rule in R language, along with its theory and
mathematical formulae.
 LAB EXERCISE 8

Aim: Monte Carlo Simulation in R

Theory

Monte Carlo simulation is a method used for sensitivity analysis, examining how a model
reacts to randomly generated inputs. The process involves three main steps: generating "N"
random inputs (often termed scenarios), executing a simulation for each input using a
computerized model, and then consolidating and evaluating the simulation outputs. Common
assessments include the mean output value, the distribution of outputs, and extreme output
values. This technique is applicable to various systems such as financial, physical, and
mathematical models. Monte Carlo simulations are independent, allowing for efficient use of
parallel computing methods, which significantly reduces computation time.

Code

num_dots = 1000

num_dots_in_circle = 0

x_inside_circle = numeric(num_dots)

y_inside_circle = numeric(num_dots)

for (i in 1:num_dots) {

x = runif(n = 1, min = -1, max = 1)

y = runif(n = 1, min = -1, max = 1)

if (x^2 + y^2 <= 1) {

num_dots_in_circle <- num_dots_in_circle + 1

x_inside_circle[num_dots_in_circle] = x

y_inside_circle[num_dots_in_circle] = y

print(4 * (num_dots_in_circle / num_dots))

plot(x_inside_circle[1:i], y_inside_circle[1:i], xlim = c(-1, 1), ylim = c(-1, 1))

[Link](0.5)

Output

Conclusion

Implemented Monte Carlo Simulation in R language, along with its theory and mathematical
formulae.
 LAB EXERCISE 9

Aim: Random Number Generation in R

Theory

At the hearth of any simulation model, there is the capability of creating numbers that mimic
those we would expect in real life. In simulation modeling we will assume that specific
processes will be distributed according to a specific random variable. For instance we will
assume that an employee in a donut shop takes a random time to serve customers distributed
according to a Normal random variable with mean μ and variance σ2. In order to then carry
out a simulation the computer will need to generate random serving times. This corresponds
to simulating numbers that are distributed according to a specific distribution.

Code

x= runif(n=10,min=1,max=5)

print(x)

#Normally Distributed Random Numbers

z=rnorm(n=5,mean=0,sd=1)

print(z)

#Binomial Random Numbers

w=rbinom(n=10,size=100,prob = 0.7)

print(w)

#Sample function

rn = sample(5:20, 5)

print(rn)
Output

Conclusion

Implemented Random Number Generation in R language.

Aim: Population Growth Model in R

Theory

Population growth models are used to model population dynamics. The model is made by
deciding if the population has an exponential growth rate or a logistic growth rate based on
the nature of the environment the population grows. From there, the model is made by
plugging in known values to solve for unknowns. Some models represent growth without
environmental constraints, while others include "ceilings" determined by limited resources.

Code

initial=[Link](readline(prompt = "Input initial population: "))

n=[Link](readline(prompt = "Input Years: "))

rate=[Link](readline(prompt = "Input rate of increment/growth: "))

f=function(initial,n,rate){

ps=0

population=c(1:(n+1))

for(i in 1:n){

population[i]=initial

ps=initial*(1+rate/100)

initial=ps;

population[n+1]=initial

plot(0:n,population,main="Population
Model",xlab="years",ylab="Population",col="blue",type="l")

return(ps)

print(f(initial,n,rate)
Output

Conclusion

Implemented Population Growth Model in R language.

Aim: Prey-Predator Model in R

Theory

Modeling and Analysis of Predator-Prey Model with Fear Effect in Prey and Hunting
Cooperation among Predators and Harvesting.

Code

f=function(iTiger,iDeer,alpha,beta,gamma,delta,time){

resTiger=0

resDeer=0

popTiger=c(1:(time+1))

popDeer=c(1:(time+1))

for(i in 1:time){

popTiger[i]=iTiger

resTiger=iTiger*(1-alpha/100)+(beta/100)*iDeer

iTiger=resTiger

popDeer[i]=iDeer

resDeer=-1*iTiger*(delta/100)+(1+gamma/100)*iDeer

iDeer=resDeer

popTiger[time+1]=iTiger

popDeer[time+1]=iDeer

plot(0:time,popTiger,main="Prey Predator
Model",xlab="years",ylab="Population",col="red",type="l")

lines(0:time,popDeer,col="black")
 legend(500,450,legend=c("Tiger Pop","Deer Pop"), col=c("red","black"),

lty=c(1,2), ncol=1)

f(500,200,0.5,0.4,0.1,0.17,1000)

Output

Conclusion

Implemented Prey-Predator Model in R language.

Aim: Implement Queuing Model in R

Theory

There is a whole area of probability called, queuing theory, which studies the mathematical
foundations and properties of such models. Queuing Theory stands as a fundamental pillar in
the domain of Modeling and Simulation, specifically focusing on systems where entities -
often customers or tasks - wait in lines for service. This mathematical discipline delves into
analyzing and predicting the behavior of queues, encompassing a wide array of applications
from telecommunications and transportation systems to healthcare and manufacturing. By
examining essential parameters like arrival rates, service times, and the number of servers,
Queuing Theory enables the quantification and understanding of system performance, aiding
in the optimization of resource allocation, reducing wait times, and enhancing overall system
efficiency. Its application in Modeling and Simulation offers a powerful tool for simulating,
evaluating, and fine-tuning complex systems, contributing to better-informed
decision-making and design across diverse industries.

Code

queue<-function(lam,mhu,n){

rho=lam/mhu

p0=1-rho

pn=(rho^n)*(1-rho)

customers=c(0)

probability=c(p0)

for(i in 1:n){

pi=(rho^i)*(1-rho)

customers=append(customers,i)

probability=append(probability,pi)

L=rho/(1-rho)

Lq=(rho^2)/(1-rho)
 W=L/lam

Wq=Lq/lam

print(c(rho,p0,pn,L,Lq,W,Wq))

plot(customers,probability,type="l",col="black",xlab="Customers",ylab="Probability")

queue(30,40,200)

Output

Conclusion

Implemented Queuing Model in R language.