EAB 2113 NUMERICAL METHOD
Question 1
Develop an M-file to implement the bisection method. Using this program solve the following problem. The velocity of falling parachutist is given as .
v (t ) = Where
gm c
(1 − e −
( c / m ) t
= velocity of parachutist parachutist =
v (t )
,
40 m / s
= gravitational constant =
g
)
,
9.8 m / s 2 = the mass of the parachutist =
m
.
68.1 kg
% DataFind obtain from coefficient, the question the drag c at the time
seconds using the initial bracket
t = 10 f=inline('((9.81*68.1)/x)*(1-exp(-(x/68.1)*10))-40' f=inline( '((9.81*68.1)/x)*(1-exp(-(x/68.1)*10))-40' ,'x' 'x'); ); xl=13; %. xu=16;of the root as [13, 16] and iterate until es=0.001; 00 1 ε ≤ 0.001 a
%computation xr=xl; %initiation while(1) while (1) %since we dont know how many interation will take place xrold=xr; %keep the previous xr xr=(xl+xu)/2; %formula for bisection method if xr~=0, ea=abs((xr-xrold)/xr)*100; end %calculate ea if xr not real answer test=f(xl)*f(xr); if test<0 xu=xr; elseif test>0 xl=xr; else ea=0; end if ea<=es, break break, ,end %end loop if the situation is satisfied end format long long; ; disp('the disp('the root for this equation is : ' ) disp(xr)
1
EAB 2113 NUMERICAL METHOD
2
EAB 2113 NUMERICAL METHOD
Question 2
Develo Develop p an M-file M-file to implem implement ent the false false posit position ion method method.. solve solve the follow following ing problem. The velocity of falling parachutist is given as .
v (t ) =
gm c
(1 − e −
( c / m ) t
Where
)
= velocity of parachutist parachutist =
v (t )
,
40 m / s
= gravitational constant =
g
,
9.8 m / s 2
3
EAB 2113 NUMERICAL METHOD
= the mass of the parachutist =
m Find the drag coefficient,
.
68.1 kg c
at the time
seconds using the initial bracket
t = 10 of the root as [13, 16] and iterate until
%. ε a
00 1 ≤ 0.001
f=inline('((9.81*68.1)/x)*(1-exp f=inline( '((9.81*68.1)/x)*(1-exp (-(x/68.1)*10))-40' ,'x' 'x'); ); xl=13; xu=16; es=0.001; xr=xl while (1) xrold=xr; xr=xu-(((f(xu)*(xl-xu))/f(xl)-f(xu))); if xr~=0, ea=abs((xr-xrold)/xr)*100; end test=f(xl)*f(xr); if test<0 xu=xr; elseif test>0 xl=xr; else test>0 ea=0; end if ea<=es, break break, , end end format long long; ; disp ('the ('the root for this equation is : ' ) disp(xr)
4
EAB 2113 NUMERICAL METHOD
5
EAB 2113 NUMERICAL METHOD
Question 3
Locate the root of
f ( x ) = 2 sin( x ) − x (a)
Using th the MA MATLAB fu function fzero with an initial guess of
.
x 0 = 2 (b)
Using Using Newton Newton-Ra -Raphs phson on method method by writ writing ing a func functio tion n M-fi M-file. le. Use an initial guess of
and iterate until
x 0 = 0.5
%. ε a
≤ 0.001
a)
6
EAB 2113 NUMERICAL METHOD
function root=newraph1(f,df,xr,es) %f=inline('2*sin(sqrt(x))-x') %df=inline('-cos(sqrt(x))/sqrt(x)-1') while(1) while (1) xr_old=xr; xr=xr-(f(xr)/df(xr)); if xr~=0,ea=abs((xr-xr_old)/xr)*100; end if ea<=es,break ea<=es, break, ,end end fprintf('\n\nthe fprintf( '\n\nthe root for this question is %2.6f\n' ,xr);
7
EAB 2113 NUMERICAL METHOD
Question 4
Develop an M-file to implement the modified secant method. Using this program determine determine the loest positive root of
with an
f ( x) = 8 sin( x) e − x f=inline('8*sin(x)*exp(-x)-1' f=inline( '8*sin(x)*exp(-x)-1' ,'x' 'x') ) Ea=1; initial guess of and . Iterate until xo=0.3; δ = 0.01 x0 = 0.3 is %.2f' ε fprintf('The fprintf( 'The initial guess ,xo) a
−1 .
= 0.000001%
while(Ea>0.000001) while(Ea>0.000001) xold=xo; xo=xo-(0.01*xo*f(xo))/(f(xo+0.01+xo)-f(xo)); Ea=((abs(xo-xold))/xo)*100; end fprintf('\nThe fprintf( '\nThe lowest positive root of function is %10.6f with approximate error %7.4f',xo,Ea) %7.4f' ,xo,Ea)
8
EAB 2113 NUMERICAL METHOD
Question 5
Find the solution of the following set of linear algebraic equations
9
EAB 2113 NUMERICAL METHOD
x + 2 y + 3 z = 1 3 x + 3 y + 4 z = 1 2 x + 3 y + 3 z = 2
a) Using the left division
b) Using Gaussian elimination
10
EAB 2113 NUMERICAL METHOD
c) Using LU decomposition
11
EAB 2113 NUMERICAL METHOD
Question 6
Develop a function M-file Tridiag.m to solve the following tridiagonal system with the Thomas algorithm. f 1 g 1 e f g 2 2 2 e3 f 3 g 3 . . .
. .
.
.
.
.
en −1
f n −1 en
x1 r 1 x r 2 2 x 3 r 3 × . = . . . . . g n −1 x n −1 r n −1 f n x n r n
Thomas Algorithm: (i)
Decomposition: 12
EAB 2113 NUMERICAL METHOD
and
ek =
ek f k −1 , where
f k (ii)
k = 2,3,4,− − − − −− , n
= f k − ek . g k −1
Forward substitution: , where
r k (iii)
.
.
k = 2,3,4,− − − − −− , n
= r k − ek . r k −1
Back substitution:
x n =
r n f n
and
, where
x k =
(r k − g k . x k +1 )
.
k = n − 1, n − 2,− − −− ,2,1
f k
Using your program, solve the following tridiagonal system.
2.01475 − 0.020875 × x1 − 0.020875 2.01475 − 0.020875 x2 − 0.020875 2.01475 − 0.020875 x3 function x=Tri(e,f,g,r) − 0.020875 2.01475 x4 %e=input('e= ');
=
4.175 0 0 2.0875
%f=input('f= '); %g=input('g= '); %r=input('r= '); e=[0 -0.020875 -0.020875 -0.020875]; f=[2.01475 2.01475 2.01475 2.01475]; g=[-0.020875 -0.020875 -0.020875 0]; r=[4.175 0 0 2.0875]; n=length(f); for k=2:n factor=e(k)/f(k-1); f(k)=f(k)-factor*g(k-1); r(k)=r(k)-factor*r(k-1); fprintf('factor=%2.4f\t fprintf( 'factor=%2.4f\t f(k)=%2.4f\t end x(n)=r(n)/f(n); for k=n-1:-1:1 x(k)=(r(k)-g(k)*x(k+1))/f(k); fprintf('x(k)=%2.4f\n' fprintf( 'x(k)=%2.4f\n',x(k)) ,x(k)) end
r(k)=%2.4f\n' ,factor,f(k),r(k))
13
EAB 2113 NUMERICAL METHOD
14
EAB 2113 NUMERICAL METHOD
Question 7
Q7.
Deve Develo lop p a MATL MATLAB AB scrip scriptt file file to deter determi mine ne the solut solutio ion n of the the follow followin ing g syste system m of linear linear equati equations ons using using the GaussGauss-Sei Seidel del iterat iteration ion method method by perfor performin ming g first first seven seven iterations. 9 x1 − 2 x 2 + 3 x 3 + 2 x 4 = 54.5 2 x1 + 8 x 2 − 2 x3 + 3 x 4 = −14
− 3 x1 + 2 x 2 + 11 x3 − 4 x 4 = 12.5 − 2 x1 + 3 x 2 + 2 x 3 + 10 x 4 = −21
i=1;x1=0;x2=0;x3=0;x4=0; disp(' disp(' i x1 x2 x3 x4') x4' ) fprintf('%2.0f fprintf( '%2.0f %8.5f %8.5f %8.5f %8.5f\n' ,i,x1,x2,x3,x4) for i=2:8 x1=(54.5-(-2*x2+3*x3+2*x4))/9; x2=(-14-(2*x1-2*x3+3*x4))/8; x3=(12.5-(-3*x1+2*x2-4*x4))/11; x4=(-21-(-2*x1+3*x2+2*x3))/10; fprintf('%2.0f fprintf( '%2.0f %8.5f %8.5f %8.5f %8.5f\n' ,i,x1,x2,x3,x4) end
15
EAB 2113 NUMERICAL METHOD
Question 8
.
Develop a script M-file to estimate
using Lagrange interpolating polynomials of f (2.75)
order 1, 2 and 3 for the following data. x f ( x )
0 0
1 0 .5
2 0 .8
3 0 .9
4 0.941176
5 0.961538
For each estimate find the true percent relative error if the try function is given by
.
f ( x ) =
x 2 (1 + x 2 )
function fint=LagrangeINT(x,y,xint) n=length(x); for i=1:n L(i)=1; for j=1:n if j~=i L(i)=L(i)*(xint-x(j))/(x(i)-x(j)); end end end f=(xint^2)/(1+xint^2); et=abs((f-sum(y.*L))/f)/100; fprintf('\nx=%.8f\n' fprintf( '\nx=%.8f\n',sum(y.*L)); ,sum(y.*L)); fprintf('Error=%6f%%\n' fprintf( 'Error=%6f%%\n',et); ,et);
16
EAB 2113 NUMERICAL METHOD
Question 9
.
The The forc force e on on a sai sailb lboa oatt mas mastt can can be repr repres esen ente ted d by by the the foll follow owin ing g func functi tion on::
F = where
z =
Compute
H
∫ 0
z e −2.5 z / H dz 7 + z
200 20 0
the the elevat elevation ion above above the deck deck and
F
the the heigh heightt of the the mast. mast.
H = for the case where
using
H = 30 17
EAB 2113 NUMERICAL METHOD
(i) (i)
the the M-f M-fil ile e for for Trap Trapez ezoi oida dall rul rule e wit with h the the step step size ize
.
h = 0.1
the MATLAB trapz function.
function evaluate = trapezoidal1(f,H,h) %f=inline('200*(z/(7+z))*exp(-25*z/30)') %y(for the trapz function)=200*(z./(7+z)).*exp(-25.*z./30) a=0; while(1) while (1) n=(H-a)/h; t0=f(a); t1=0; for i=1:n-1 t1=t1+f(a+h*i); end t2=f(H); value=h/2*(t0+2*t1+t2); break end fprintf('The fprintf( 'The Trapezoidal Integral: %2.8f\n' ,value)
Question 10
18
EAB 2113 NUMERICAL METHOD
Develop an M-file to implement Simpson’s 1/3 rule. Using your program solve the following problem. The velocity of falling parachutist is given as .
v(t ) = Where
gm c
(1 − e −( c / m ) t )
= velocity of parachutist, parachutist,
v (t ) = gravitational constant =
g
,
9.8 m / s 2 = mass of the parachutist parachutist =
m
,
45 kg = the drag coefficient =
c
.
32.5 kg / s
If the distance, d, traveled by the parachutist is given by , 6
∫
d = v(t )dt 0
find the distance using Simpson’s 1/3 rule for the segments 10, 20, 50, and 100.
function distance=simpson(V,a,b,n) h=(b-a)/n; t(1)=a; for i=2:n+1 t(i)=t(i-1)+h; end P=0; for i=2:2:n P=P+V(t(i)); end T=0; for i=3:2:n-1 T=T+V(t(i)); end d=(h/3)*(V(a)+(4*P)+(2*T)+V(b)); fprintf('\nThe fprintf( '\nThe distance travelled by the parachutist=%9.6f\n' ,d)
19
EAB 2113 NUMERICAL METHOD
20
EAB 2113 NUMERICAL METHOD
Question 11
Develop an M-file for Euler’s method to solve a first order ordinary differential differential equation (ODE). E i(t )
T The he curr curren entt arou around nd the the circ circui uitt at time time
is govern governed ed by the follow following ing
t differential equation ,
3
di dt
= 2i + 3e −2t
.
i (0 ) = 2
Using your program, solve the above initial value problem over the interval from
to 2 with the step size
t = 0
.
h = 0.1 21
EAB 2113 NUMERICAL METHOD
function [t,i]= Eulode(didt,tspan,i0,h) %input: %didt= name of the function f(t,i) %tspan= [ti,tf] where ti and tf= initial and final valus of independent %variable %y0= initial value of the dependent variable %h=step size %output: %[t,i] where t= vector of the independent variable % i= vector of the solution for the dependent variable tx=tspan(1); ty=tspan(2); t=(tx:h:ty)'; n=length(t); i=i0*ones(n,1); for x=1:n-1 i(x+1)=i(x)+feval(didt,t(x),i(x))*h; end
22
EAB 2113 NUMERICAL METHOD
23
EAB 2113 NUMERICAL METHOD
Question 12
Develop an M-file for Fourth-Order Runge-Kutta method to solve a first order ordinary differential equation (ODE). Using your program solve the following initial value problem over the interval from
to 2 with the step size
x = 0
.
h = 0.2 .
dy dx
=
x 2
+ y ,
y (0) = 0.8
function [x,y]=rk40de(dydx,xspan,y0,h) xi=xspan(1); xf=xspan(2); x=(xi:h:xf)'; n=length(x); y=y0*ones(n,1); for i=1:n-1 k1=feval(dydx,x(i),y(i)); k2=feval(dydx,x(i)+(1/2)*h,y(i)+(1/2)*k1*h); k3=feval(dydx,x(i)+(1/2)*h,y(i)+(1/2)*k2*h); k4=feval(dydx,x(i)+h,y(i)+k3*h); y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4)*h; end
24
EAB 2113 NUMERICAL METHOD
25