Find a real root of \(f(x)=x^3-x-1\) in the interval (1,2) by bisection method.
## [1] -1 5 -5Regula Falsi Method
Iteration Method
Find a real root of \(f(x)=x^3-x-1\) near \(x=2\) by iteration method.
Sol: Rewrite the equation as \(x=(x+1)^{1/3}\) so that \(\phi(x)=(x+1)^{1/3}\).
Newton Raphson Method
\[ x_{n+1}=x_n-\frac{f(x_n)}{f^\prime(x_n)} \]
Find a real root of \(f(x)=x^3-x-1\) near \(x=2\).
Sol: Differentiating \(f^\prime(x)=3x^2-1\) and thus the NR formula becomes \[ \begin{array}{rcl} x&=& x-\frac{x^3-x-1}{3x^2-1}\\ &=& \frac{3x^3-x-(x^3-x-1)}{3x^2-1}\\ &=& \frac{2x^3+1}{3x^2-1} \end{array} \]
Find a real root of \(f(x)=x e^x-1\) in the interval (0,1) by bisection method.
## [1] -1.000000 1.718282 -1.718282Regula Falsi Method
Iteration Method
Find a real root of \(f(x)=x e^x-1\) near \(x=0\) by iteration method.
Sol: Rewrite the equation as \(x=e^{-x}\) so that \(\phi(x)=e^{-x}\).
Newton Raphson Method
\[ x_{n+1}=x_n-\frac{f(x_n)}{f^\prime(x_n)} \]
Find a real root of \(f(x)=x e^x -1\) near \(x=0\).
Sol: Differentiating \(f^\prime(x)=e^x(x+1)\)
Find a real root of \(f(x)=cos(x)-3x\) in the interval (0,1) by bisection method.
## [1] 1.000000 -2.459698 -2.459698Regula Falsi Method
Iteration Method
Find a real root of \(f(x)=cos(x)-3x\) near \(x=0\) by iteration method.
Sol: Rewrite the equation as \(x=e^{-x}\) so that \(\phi(x)=e^{-x}\).
Newton Raphson Method
\[ x_{n+1}=x_n-\frac{f(x_n)}{f^\prime(x_n)} \]
Find a real root of \(f(x)=cos(x)-3x\) near \(x=0\).
Sol: Differentiating \(f^\prime(x)=-sin(x)-3\)
Newton Raphson Method from Taylor’s Series:
Expanding the given function \(f(x)\) by Taylor’s series about \(x=a\), \[ f(x)=f(a)-(x-a)f^\prime(a)+O(h^2) \] Let x be such that \(f(x)=0\) so that \[ x=a-\frac{1}{f^\prime(a)} f(a) \] Setting \(x=x_1\) and \(a=x_0\), we get the Newton Raphson formula.
Newton Raphson Method to Solve simultaneous equations:
Find one real root of the following equations \[ \begin{array}{c} f(x,y)=0\\ g(x,y)=0 \end{array} \] near \(X_0=(x_0,y_0)\)
Expanding the given functions in Taylor Series about \(X_0\), \[ \begin{array}{c} f(x,y)=f(x_0,y_0)+(x-x_0)f_x(x_0,y_0)+(y-y_0)f_y(x_0,y_0)\\ g(x,y)=g(x_0,y_0)+(x-x_0)g_x(x_0,y_0)+(y-y_0)g_y(x_0,y_0) \end{array} \] the next approximate root will be \[ \begin{array}{lcr} \left[ \begin{array}{c} x_1\\y_1 \end{array} \right] &=& \left[ \begin{array}{c} x_0\\y_0 \end{array} \right]-\left[ \begin{array}{cc} f_x^0 & f_y^0\\g_x^0&g_y^0 \end{array} \right]^{-1}\left[ \begin{array}{c} x_0\\y_0 \end{array} \right] \end{array} \] So the Jacobian matrix (J) is to be computed each time and inverted to get the next approximated root.
Find the solution of the simultaneous equations \[ \begin{array}{c} x^2+xy+y^2=7\\ x^3+y^3=9 \end{array} \] near \(x_0=1.5\) and \(y_0=0.5\).
Ans: Let \(f(x,y)=x^2+xy+y^2-7=0\) and \(g(x,y)=x^3+y^3-9=0\). Differentiating partially w.r.t. x and y we get \[ \begin{array}{l} f_x=2x+y\\ f_y=x+2y\\ g_x=3x^2\\ g_y=3y^2 \end{array} \] The Jacobian of the functions is \[ J\left(\frac{f,g}{x,y} \right) \] The inverse of the Jacobian is \[ J^{-1}=\frac{1}{f_x^0g_y^0-g_x^0f_y^0}\left[ \begin{array}{rr} g_y^0 & -f_y^0\\ -g_x^0 & f_x^0 \end{array} \right] \] It is to be calculated for each iteration and next approximation is given by \[ X_1=X_0-J^{-1}X_0 \]
## [1] "the Jacobian is"## [,1] [,2]
## [1,] 3.50 2.50
## [2,] 6.75 0.75## [1] "inverse Jacobian is"## [,1] [,2]
## [1,] -0.05263158 0.1754386
## [2,] 0.47368421 -0.2456140## [1] "the next approximation root is"## [,1]
## [1,] 2.2675439
## [2,] 0.9254386## [,1]
## [1,] 2.2675439
## [2,] 0.9254386
## [,1]
## [1,] 2.0372712
## [2,] 0.9644695
## [,1]
## [1,] 2.0012578
## [2,] 0.9987366Power method
Power method is used to find the dominant eigen vector and hence it’s eigen value. The initial vector is multiplied by the matrix A successively simultaneously dividing the vector by the largest value of the vector. The procedure is continued two consecutive vectors do not differ significantly.
Find the dominant eigen vector and hence eigen value of the following matrix \[\left[ \begin{array}{rrr} 1 & 2 & 0\\ -2 & 1 & 2\\ 1 & 3 & 1 \end{array}\right] \] starting with the initial vector x0 = (1,1,1).
x0 -> a x0 -> x1=a x0/max(a x0) -> x0=x1 -> …. -> xn
## [1] "initial vector x0"## x0
## 1 1
## 2 1
## 3 1## [1] "the vector a x0"## a.....x0
## 1 3
## 2 1
## 3 5## [1] "the next vector ax0*x0/max(ax0*x0)"## a.....x0
## 1 0.6
## 2 0.2
## 3 1.0## [1] "dominant eigen vector is"## [1] "max eigen value is ax0*x0/(x0*x0)" "3"Jacobi Method Solve the equations \(2x+y+6z=9\), \(8x+3y+2z=13\), \(x+5y+z=7\)
Using Jacobi method Step 1: Re arrange the equations \(8x+3y+2z=13\), \(x+5y+z=7\), \(2x+y+6z=9\) Step 2: Use 1st equation for 1st variable, 2nd for 2nd and so on x= 1/8 (13-3y-2z) y=1/5(7-x-z) z=1/6 (9-2x-y) Step 3: Let x0,y0,z0 initial approximation the next approximation will be x1= 1/8 (13-3y0-2z0) y1=1/5(7-x0-z0) z1=1/6 (9-2x0-y0) Step 4: if X0-X1 is small stop the solution is x1,y1,z1 Step 5: x0=x1, y0=y1 and z0=z1
## x y z
## [1,] 1.625 1.400 1.500
## [2,] 0.725 0.775 0.725
## [3,] 1.153 1.110 1.129
## [4,] 0.926 0.944 0.931
## [5,] 1.039 1.029 1.034
## [6,] 0.981 0.986 0.982
## [7,] 1.010 1.007 1.009
## [8,] 0.995 0.996 0.995
## [9,] 1.003 1.002 1.002
## [10,] 0.999 0.999 0.999Gauss Seidal Method Step 1: Re arrange the equations \(8x+3y+2z=13\), \(x+5y+z=7\), \(2x+y+6z=9\) Step 2: Use 1st equation for 1st variable, 2nd for 2nd and so on x= 1/8 (13-3y-2z) y=1/5(7-x-z) z=1/6 (9-2x-y) Step 3: Let x0,y0,z0 initial approximation the next approximation will be x1= 1/8 (13-3y0-2z0) y1=1/5(7-x1-z0) z1=1/6 (9-2x1-y1) Step 4: if X0-X1 is small stop the solution is x1,y1,z1 Step 5: x0=x1, y0=y1 and z0=z1
## x y z
## [1,] 1.625 1.075 0.779
## [2,] 1.027 1.039 0.985
## [3,] 0.989 1.005 1.003
## [4,] 0.997 1.000 1.001Using Gauss Seidel method solve the equations \(10x+y+z=12\) , \(2x+10y+z=13\), \(2x+2y+10z=14\)
x0=0
y0=0
z0=0
df=c()
for(i in 1:10)
{
x1=(12-y0-z0)/10
y1=(13-2*x1-z0)/10
z1=(14-2*x1-2*y1)/10
df=rbind(df,c(x1,y1,z1))
if(max(abs(c(x1,y1,z1)-c(x0,y0,z0)))<=0.01) break
x0=x1
y0=y1
z0=z1
}
colnames(df)<-c('x','y','z')
round(df,3)## x y z
## [1,] 1.200 1.060 0.948
## [2,] 0.999 1.005 0.999
## [3,] 1.000 1.000 1.000The problem is to evaluate the definite integrals \[ \int_a^bf(x)\, dx \]
Here one of the following are to be specified \begin{enumerate} the step size (h) no of division points (n)ready made table is specified (xi,fi), i=1,n \end{}enumerate)
Case-1
If step size is given find the division points by adding h to a and successively there after
a=0
b=1
h=0.2
print("points of division")## [1] "points of division"seq(a,b,h)## [1] 0.0 0.2 0.4 0.6 0.8 1.0Case-2 If the no of divisions are given then calculate h=(b-a)/n and follow case-1.
a=0
b=1
n=4
h=(b-a)/n
print("points of division")## [1] "points of division"seq(a,b,h)## [1] 0.00 0.25 0.50 0.75 1.00Case-3 A tabulated function is given as follows
f=function(x) exp(-x)
x=seq(0,1,0.2)
data.frame(x,round(f(x),3))\[ \int_{x_0}^{x_n} f(x)dx=\frac{h}{2}\left[f_0+2(f_1+f_2+\cdots+f_{n-1})+f_n \right] \] here \(a=x_0\) and \(b=x_n\) taking \(n\) number of divisions.
## [1] "the value of the integral is" "0.6342"#input data initial, final and no of intervals
a=0
b=1
n=6 #enter any even no
#prepare the points of division
h=(b-a)/n
x=seq(0,1,h)
x## [1] 0.0000000 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333 1.0000000#define the function
fx=round(exp(-x),3)
#prepare the weights
w=c()
for (i in 2:n-1)
{
temp=if(i%%2==1) 4 else 2
w=append(w,temp)
}
w=c(1,w,1)
w## [1] 1 4 2 4 2 4 1#prepare the table
#w=c(1,rep(2,length(x)-2),1)
df=data.frame(x,fx,w,w*fx)
df#use the formulae
res=h/3*sum(df$w...fx)
c('the value of the integral is',round(res,4))## [1] "the value of the integral is" "0.6322"\[ \int_{x_0}^{x_1}f(x)\, dx=\frac{3h}{8} \left[x_0+x_n+3(y_1+y_2+\cdots)+2(y_3+y_6+\cdots multiples\, of\, 3 ) \right] \]
## [1] 0.00 0.17 0.33 0.50 0.67 0.83 1.00## [1] 1 3 3 2 3 3 1## [1] "the value of the integral is" "0.6322"\[ y^\prime=x^2-y,\: y(0)=1,\: h=0.1 \] Find the value of \(y(0.1)\) correct up to four decimal places.
f=function(x,y) x^2-y
x0=0
y0=1
h=0.1
y10=y0+h*f(x0,y0)
#y10
for (i in 1:3)
{
y11=y0+h/2*(f(x0,y0)+f(x1,y10))
print(round(y11,4))
y10=y11
x1=x0+h
}## [1] 0.955
## [1] 0.9028
## [1] 0.9054## [1] 4 3 4 6 8
## [1] -1 1 2 2
## [1] 2 1 0
## [1] -1 -1
## [1] 0## [1] "final"## [1] 6## [1] 4 -1 2 -1 0## [1] 8 2 0 -1 0\[ y^\prime =y-x,\, y(0)=2,h=0.2 \] Find \(y(0.2)\).
f=function(x,y) y-x
x0=0
y0=2
h=0.2
k1=h*f(x0,y0)
k2=h*f(x0+h/2,y0+k1/2)
k3=h*f(x0+h/2,y0+k2/2)
k4=h*f(x0+h,y0+k3)
res=y0+1/6*(k1+2*k2+2*k3+k4)
data.frame(k1,k2,k3,k4,res)\[ \begin{array}{|c|c|c|c|c|c|c|} \hline x&y &dy & d^2y &d^3y &d^4y &d^5y\\ \hline 0&7&&&&& \\ &&4&&&& \\ 5&11&&-1 &&&\\ & &3& &2&&\\ 10 &14 & &1& &-1&\\ & &4 & &1 & &0\\ 15&18 & &2 & &-1 &\\ & &6 & &0 & &\\ 20 &24 & &2 & & &\\ & &8 & & & &\\ 25&32 & & & & &\\ \hline \end{array} \]
xcx