OUTPUT: (1)
5
-1
6
0.6666667
0.6931472
0.1353353
-1
-2.449213e-16
2
-1+0i
-2-2i
-5+10i
0.44+0.08i
0.5+0.5i
3
2
3.605551
0.5880026
3-2i
0.7071068+0.7071068i
0.7853982
-11.3
4
-11.3
-0.5116667
5.580288
31.13962
6
9.60 5.10 -33.90 -2.01 12.00 0.00
-0.02919952 0.65998315 0.80614681 0.94441031 -0.41614684 1.00000000
Min. 1st Qu. Median Mean 3rd Qu. Max.
-11.3000 -0.5025 0.8500 -0.5117 2.8250 4.0000
0
[,1] [,2]
[1,] 4 8
[2,] 12 16
[,1] [,2]
[1,] -2 -4
[2,] -6 -8
[,1] [,2]
[1,] 21 30
[2,] 45 66
[,1] [,2]
[1,] 1 3
[2,] 2 4
-2
OUTPUT: (3)
"Gauss seidel method converged in 5 iterations"
3.0166004 1.9859643 0.9118754
BISECTION METHOD GAUSS SEIDEL METHOD
bisectionroot=function(f,xmin,xmax,tol=1e-5) gauss_seidel=function(A,b,x0,tol=1e-
{ 3,max_iter=100)
a=xmin; b=xmax {
if(a>=b) n=length(b)
{ x=x0
cat("error:xmin>xmax\n") for(iter in 1:max_iter)
return(NULL) {
} x_old=x
if(f(a)==0) for(i in 1:n)
{ {
return(a) sum_Ax=0
} for(j in 1:n)
else if(f(b)==0) {
{ if(i!=j)
return(b) {
} sum_Ax=sum_Ax+A[i,j]*x[j]
else if(f(a)*f(b)>0) }
{ }
cat("error:f(xmin) and f(xmax) have the x[i]=(b[i]-sum_Ax)/A[i,i]
same sign\n") }
return(NULL) if(max(abs(x-x_old))<tol)
} {
iter=0 print(paste("Gauss seidel method
while((b-a)>tol) converged in",iter,"iterations"))
{ return(x)
c=(a+b)/2 }
if(f(c)==0) }
{ print("Gauss seidel method did not converge
return(c) within the maximum iterations")
} return(x)
else if(f(a)*f(c)<0) }
{ A=matrix(c(8,-3,2,4,11,-1,6,3,12),
b=c nrow=3,byrow=TRUE)
} b=c(20,33,35)
else x0=c(0,0,0)
{ gauss_seidel(A,b,x0)
a=c
} OUTPUT: (2)
iter=iter+1 8.028069e-01 1.600000e+01 7.629395e-06
}
return(c((a+b)/2,iter,(b-a)))
}
f=function(x)x^3-sin(x)^2
bisectionroot(f,0.5,1)
INTRODUCTION
2+3 CRAMER’S RULE
2-3 determinant_matrix=function(A)
2*3 {
2/3 return(det(A))
log(2) }
exp(-2) cramers_rule=function(A,b)
cos(pi) {
sin(2*pi) n=nrow(A)
abs(-2) det_A=determinant_matrix(A)
(1i)^2 if(det_A==0)
(1+2i)-(3+4i) {
(1+2i)*(3+4i) stop("The determinant of A is zero,the system
(1+2i)/(3+4i) has no unique solution")
(1i+(1i)^2+(1i)^3+(1i)^4+(1i)^5)/(1+1i) }
Re(3+2i) solution=numeric(n)
Im(3+2i) for(i in 1:n)
Mod(3+2i) {
Arg(3+2i) A_i=A
Conj(3+2i) A_i[,i]=b
exp(pi*1i/4) solution[i]=(determinant_matrix(A_i))/det_A
x=pi/4 }
y=sin(x) return(solution)
asin(y) }
x=c(3.2,1.7,-11.3,-0.67,4,0) A=matrix(c(10,-5,-2,4,-
x[3] 10,3,1,6,10),nrow=3,byrow=TRUE)
max(x) b=c(3,-3,-3)
min(x) determinant_matrix(A)
mean(x) cramers_rule(A,b)
sd(x)
var(x) SIMPSON’S 3/8 RULE
length(x) f=function(x)
3*x {
cos(x/2) 1/(1+(x*x))
summary(x) }
prod(x) a=0
m=matrix(c(1,2,3,4), nrow=2, byrow=TRUE) b=6
n=matrix(c(3,6,9,12), nrow=2, byrow=TRUE) n=6
m+n h=(b-a)/n
m-n sum=0
m%*%n for(i in 1:(n-1))
t(m) {
det(m) x=a+(i*h)
if(i%%3==0)
{
OUTPUT: (4) sum=sum+2*f(x)
-1063 }
0.3414864 0.2850423 -0.5051740 else
{
sum=sum+(3*f(x))
}
}
value=(3*h/8)*(f(a)+f(b)+sum)
print(value)
GAUSS JACOBI METHOD TRAPEZOIDAL METHOD
jacobi=function(A,b,x0,tol=1e- f=function(x)
6,max_iter=100) {
{ x^4
n=length(b) }
x=x0 a=-3
x_new=x0 b=3
for(iter in 1:max_iter) n=6
{ h=(b-a)/n
for(i in 1:n) sum=0
{ for(i in 1:(n-1))
sum_Ax=0 {
for(j in 1:n) x=a+(i*h)
{ sum=sum+f(x)
if(i!=j) print(f(x))
{ }
sum_Ax=sum_Ax+A[i,j]*x[j] value=(h/2)*(f(a)+f(b)+2*sum)
} print(value)
}
x_new[i]=(b[i]-sum_Ax)/A[i,i] NEWTON RAPHSON METHOD
} newton_raphson=function(f,fd,x0,tol=1e-
if(max(abs(x_new-x))<tol) 5,max_iter=20)
{ {
print(paste("Jacobi method converges x=x0
in",iter,"iteration")) for(i in 1:max_iter)
return(x_new) {
} value=x-(f(x)/fd(x))
x=x_new if(abs(value-x)<tol)
cat("iteration",iter,":",x_new,"\n") {
} return(value)
print("Jacobi method does not converge in }
maximum nuber of iterations") x=value
return(x_new) }
} print("did not converge")
A=matrix(c(10,-5,-2,4, return(NA)
-10,3,1,6,10),nrow=3,byrow=TRUE) }
b=c(3,-3,-3) f=function(x)
x0=c(0,0,0) {
result=jacobi(A,b,x0) x-cos(x)
print("solution is") }
print(result) fd=function(x)
{
OUTPUT: (5) 1+sin(x)
"Jacobi method converges in 15 iteration" }
"solution is" x0=1
0.3414867 0.2850426 -0.5051745 root=newton_raphson(f,fd,x0)
print(root)
OUTPUT: (6)
16 OUTPUT: (7)
1 0.7390851
0
1 OUTPUT: (9)
16 1.357081
115
SIMPSON’S 1/3 RULE NEWTON’S FORWARD DIFFERENCE
f=function(x) FORMULA
{ x=c(3,3.2,3.4,3.6,3.8)
1/(1+(x*x)) y=c(-14,-10.032,-5.296,-0.256,6.672)
} h=x[2]-x[1]
a=0 n=length(y)
b=6 y1=numeric(n-1)
n=6 for(i in 1:(n-1))
h=(b-a)/n {
sum=0 y1[i]=y[i+1]-y[i]
for(i in 1:(n-1)) }
{ y2=numeric(n-2)
x=a+(i*h) for(i in 1:(n-2))
if(i%%2==0) {
{ y2[i]=y1[i+1]-y1[i]
sum=sum+2*f(x) }
} y3=numeric(n-3)
else for(i in 1:(n-3))
{ {
sum=sum+(4*f(x)) y3[i]=y2[i+1]-y2[i]
} }
} y4=numeric(n-4)
value=(h/3)*(f(a)+f(b)+sum) for(i in 1:(n-4))
print(value) {
y4[i]=y3[i+1]-y3[i]
OUTPUT: (8) }
1.366081 print(y)
print(y1)
RUNGE-KUTTA METHOD print(y2)
F=function(x,y) print(y3)
{ print(y4)
return(-y-(x*y*y)) dy_dx=(1/h)*(y1[1]-(y2[1]/2)+(y3[1]/3)-
} (y4[1]/4))
rk4=function(f,x0,y0,x_end,h) print(dy_dx)
{
x=x0
y=y0
while(x<x_end) OUTPUT:
{ "value of y at x= 0.1 is 0.900623706760202"
k1=h*f(x,y) "value of y at x= 0.2 is 0.804631486126234"
k2=h*f(x+0.5*h,y+0.5*k1) "value of y at x= 0.3 is 0.71443039481972"
k3=h*f(x+0.5*h,y+0.5*k2)
k4=h*f(x+h,y+k3) OUTPUT (10)
y=y+((k1+(2*k2)+(2*k3)+k4)/6)
x=x+h
print(paste("value of y at x=",x,"is",y))
}
}
x0=0
y0=1
x_end=0.3
h=0.1
result=rk4(f,x0,y0,x_end,h)
OUTPUT: (11)
-14.000 -10.032 -5.296 -0.256 6.672
3.968 4.736 5.040 6.928
0.768 0.304 1.888
-0.464 1.584
2.048
14.58667
NEWTON’S BACKWARD DIFFERENCE FORMULA
x=c(1931,1941,1951,1961,1971)
y=c(40.62,60.8,79.95,103.56,132.65)
h=x[2]-x[1]
n=length(y)
y1=numeric(n-1)
for(i in 1:(n-1))
{
y1[i]=y[i+1]-y[i]
}
y2=numeric(n-2)
for(i in 1:(n-2))
{
y2[i]=y1[i+1]-y1[i]
}
y3=numeric(n-3)
for(i in 1:(n-3))
{
y3[i]=y2[i+1]-y2[i]
}
y4=numeric(n-4)
for(i in 1:(n-4))
{
y4[i]=y3[i+1]-y3[i]
}
print(y)
print(y1)
print(y2)
print(y3)
print(y4)
dy_dx=(1/h)*(y1[n-1]+(y2[n-2]/2)+(y3[n-3]/3)+(y4[n-4]/4))
print(dy_dx)
OUTPUT: (12)
40.62 60.80 79.95 103.56 132.65
20.18 19.15 23.61 29.09
-1.03 4.46 5.48
5.49 1.02
-4.47
3.10525