Chapter 23 First Deriative

23.10 Option Explicit Sub RhombTest() Dim maxit As Integer Dim a As Single, b As Single, es As Single Dim x As Single x = 0.5 maxit = 3 es = 0.001 MsgBox RhomDiff(x, maxit, es) End Sub Function RhomDiff(x, maxit, es) Dim n As Integer, j As Integer, k As Integer, iter As Integer Dim i(10, 10) As Single, ea As Single, del As Single, a As Single, b As Single n = 1 i(1, 1) = DyDx(x, n) iter = 0 Do iter = iter + 1 n = 2 ^ iter i(iter + 1, 1) = DyDx(x, n) For k = 2 To iter + 1 j = 2 + iter - k i(j, k) = (4 ^ (k - 1) * i(j + 1, k - 1) - i(j, k - 1)) / (4 ^ (k - 1) - 1) Next k ea = Abs((i(1, iter + 1) - i(1, iter)) / i(1, iter + 1)) * 100 If (iter >= maxit Or ea <= es) Then Exit Do Loop RhomDiff = i(1, iter + 1) End Function Function DyDx(x, n) Dim a As Single, b As Single a = x - x / n b = x + x / n DyDx = (f(b) - f(a)) / (b - a) End Function Function f(x) f = -0.1 * x ^ 4 - 0.15 * x ^ 3 - 0.5 * x ^ 2 - 0.25 * x + 1.2 End Function 23.11 The following program implements Eq. 23.9. Option Explicit Sub TestDerivUnequal() Dim n As Integer, i As Integer Dim x(100) As Single, y(100) As Single, dy(100) As Single Range("a5").Select n = ActiveCell.Row Selection.End(xlDown).Select n = ActiveCell.Row - n Range("a5").Select For i = 0 To n x(i) = ActiveCell.Value ActiveCell.Offset(0, 1).Select y(i) = ActiveCell.Value ActiveCell.Offset(1, -1).Select Next i For i = 0 To n dy(i) = DerivUnequal(x(), y(), n, x(i)) Next i Range("c5").Select For i = 0 To n ActiveCell.Value = dy(i) ActiveCell.Offset(1, 0).Select Next i End Sub Function DerivUnequal(x, y, n, xx) Dim ii As Integer If xx < x(0) Or xx > x(n) Then DerivUnequal = "out of range" Else If xx < x(1) Then DerivUnequal = DyDx(xx, x(0), x(1), x(2), y(0), y(1), y(2)) ElseIf xx > x(n - 1) Then DerivUnequal = DyDx(xx, x(n - 2), x(n - 1), x(n), y(n - 2), y(n - 1), y(n)) Else For ii = 1 To n - 2 If xx >= x(ii) And xx <= x(ii + 1) Then If xx - x(ii - 1) < x(ii) - xx Then 'If the unknown is closer to the lower end of the range, 'x(ii) will be chosen as the middle point DerivUnequal = DyDx(xx, x(ii - 1), x(ii), x(ii + 1), y(ii - 1), y(ii), y(ii + 1)) Else 'Otherwise, if the unknown is closer to the upper end, 'x(ii+1) will be chosen as the middle point DerivUnequal = DyDx(xx, x(ii), x(ii + 1), x(ii + 2), y(ii), y(ii + 1), y(ii + 2)) End If Exit For End If Next ii End If End If End Function Function DyDx(x, x0, x1, x2, y0, y1, y2) DyDx = y0 * (2 * x - x1 - x2) / (x0 - x1) / (x0 - x2) _ + y1 * (2 * x - x0 - x2) / (x1 - x0) / (x1 - x2) _ + y2 * (2 * x - x0 - x1) / (x2 - x0) / (x2 - x1) End Function The result is An even more elegant approach is to put cubic splines through the data (recall Sec. 20.2 and the solution for Prob. 20.10) to evaluate the derivatives. 23.12 (a) Create the following M function: function y=f(x) y=9.8*68.1/12.5*(1-exp(-12.5/68.1*x)); Then implement the following MATLAB session: >> Q=quad('f',0,10) Q = 289.4351 (b) ∫0 (1 − e ) dt d (t ) = gm c d (t ) = gm  m − ( c / m ) t  t+ e  c  c 0 d (t ) = 9.8( 68.1)  68.1 − (12.5 / 68.1)10 68 .1  10 + e −0− = 289.4351  12.5  12.5 12.5  0 t − ( c / m )t t 10 (c) Create the following M function: >> function y=f(x) >> y=9.8*68.1/12.5*(1-exp(-12.5/68.1*x)); Then implement the following MATLAB session: >> x=[9.99 10.01] >> y=f(x) >> d=diff(y)./diff(x) d = 1.5634 (d) a (t ) = ( gm d 1 − e − (c / m)t c dt ) a(t ) = ge −( c / m )t a(t ) = 9.8e − (12.5 / 68.1)10 = −1.56337 23.13 (a) Create the following M function: function y=fn(x) y=1/sqrt(2*pi)*exp(-(x.^2)/2); Then implement the following MATLAB session: >> x=-2:.1:2; >> y=fn(x); >> Q=quad('fn',-1,1) Q = 0.6827 >> Q=quad('fn',-2,2) Q = 0.9545 Thus, about 68.3% of the area under the curve falls between –1 and 1 and about 95.45% falls between –2 and 2. (b) >> >> >> >> >> >> >> x=-2:.1:2 y=fn(x) d=diff(y)./diff(x) x=-1.95:.1:1.95 d2=diff(d)./diff(x) x=-1.9:.1:1.9 plot(x,d2,'o') Thus, inflection points (d2y/dx2 = 0) occur at –1 and 1. 23.14 (a) Create the following M function: function y=fn(x) y=1/sqrt(2*pi)*exp(-(x.^2)/2); Then implement the following MATLAB session: >> x=-2:.5:2; >> y=fn(x); >> Q=quad('fn',-1,1) Q = 0.6827 >> Q=quad('fn',-2,2) Q = 0.9545 Thus, about 68.3% of the area under the curve falls between –1 and 1 and about 95.45% falls between –2 and 2. (b) >> >> >> >> >> d=diff(y)./diff(x); x=-1.75:.5:1.75; d2=diff(d)./diff(x); x=-1.5:.5:1.5; plot(x,d2,'o') Thus, inflection points (d2y/dx2 = 0) occur at –1 and 1. 23.15 Program Integrate Use imsl Implicit None Integer::irule=1 Real::a=-1.,b=1,errabs=0.0,errrel=0.001 Real::errest,res,f External f Call QDAG(f,a,b,errabs,errrel,irule,res,errest) Print '('' Computed = '',F8.4)',res Print '('' Error estimate ='',1PE10.3)',errest End Program Function f(x) Implicit None Real:: x , f Real::pi Parameter(pi=3.1415927) f=1/sqrt(2*pi)*exp(-x**2/2) End Function Answers: x = -1 to 1: Computed = x = -2 to 2: Computed = x = -3 to 3: Computed = 0.6827 Error estimate = 4.069E-06 0.9545 Error estimate = 7.975E-06 0.9973 Error estimate = 5.944E-06 23.16 MATLAB Script: % Prob2316 Integration program a=0; b=pi/2; integral=quad('ff',a,b) end function y=ff(x); y=sin(sin(x)); >> prob2316 integral = 0.8932 23.17 MATLAB Script: %Numerical Integration of sin(t)/t = function sint(t) %Limits: a=0, b=2pi %Using the "quad" and "quadl" function for numerical integration %Plot of function t=0.01:0.01:2*pi; y=ff2(t); plot(t,y); grid %Integration format long a=0.01; b=2*pi; Iquad=quad('ff2',a,b) Iquadl=quadl('ff2',a,b) function y=ff2(t); y=sin(t)./t; MATLAB execution: >> prob2317 Iquad = 1.40815164305082 Iquadl = 1.40815163168846 23.18 %Centered Finite Difference First & Second Derivatives of Order O(dx^2) %Using diff(y) dx=0.5; y=[1.4 2.1 3.3 4.7 7.1 6.4 8.8 7.2 8.9 10.7 9.8]; dyf=diff(y); % First Derivative Centered FD using diff n=length(y); for i=1:n-2 dydxc(i)=(dyf(i+1)+dyf(i))/(2*dx); end %Second Derivative Centered FD using diff dy2dx2c=diff(dyf)/(dx*dx); fprintf('first derivative \n'); fprintf('%f\n', dydxc) fprintf('second derivative \n'); fprintf('%f\n', dy2dx2c) first derivative 1.900000 2.600000 3.800000 1.700000 1.700000 0.800000 0.100000 3.500000 0.900000 second derivative 2.000000 0.800000 4.000000 -12.400000 12.400000 -16.000000 13.200000 0.400000 -10.800000 23.19 % % % % % % Finite Difference Approximation of slope For f(x)=exp(-x)-x f'(x)=-exp(-x)-1 Centered diff. df/dx=(f(i+1)-f(i-1))/2dx + O(dx^2) Fwd. diff. df/dx=(-f(i+2)+4f(i+1)-3f(i))/2dx + O(dx^2) Bkwd. diff. df/dx=(3f(i)-4f(i-1)+f(i-2))/2dx + O(dx^2) x=2; fx=exp(-x)-x; dfdx2=-exp(-x)-1; %approximation dx=0.5:-0.01:.01; for i=1:length(dx) %x-values at i+-dx and +-2dx xp(i)=x+dx(i); x2p(i)=x+2*dx(i); xn(i)=x-dx(i); x2n(i)=x-2*dx(i); %f(x)-values at i+-dx and +-2dx fp(i)=exp(-xp(i))-xp(i); f2p(i)=exp(-x2p(i))-x2p(i); fn(i)=exp(-xn(i))-xn(i); f2n(i)=exp(-x2n(i))-x2n(i); %Finite Diff. Approximations Cdfdx(i)=(fp(i)-fn(i))/(2*dx(i)); Fdfdx(i)=(-f2p(i)+4*fp(i)-3*fx)/(2*dx(i)); Bdfdx(i)=(3*fx-4*fn(i)+f2n(i))/(2*dx(i)); end dx0=0; plot(dx,Fdfdx,'--',dx,Bdfdx,'-.',dx,Cdfdx,'-',dx0,dfdx2,'*') grid title('Forward, Backward, and Centered Finite Difference approximation - 2nd Order Correct') xlabel('Delta x') ylabel('df/dx') gtext('Centered'); gtext('Forward'); gtext('Backward') F orward, B ac k ward, and Centered F inite Differenc e approx im ation - 2nd O rder Correc t -1.11 -1.115 -1.12 -1.125 df/dx B ac k ward -1.13 F orward -1.135 -1.14 Centered -1.145 0 0.05 0.1 0.15 0.2 0.25 Delta x 0.3 0.35 0.4 0.45 0.5 23.20 23.21 a) v= a= dx x (ti +1) − x( ti −1) 7 .3 − 5.1 = x ′(ti ) = = = 1.1 m / s dt 2h 2 d 2x dt 2 = x′′(t i ) = x(ti +1 ) − 2 x (ti ) + x(ti −1 ) h 2 = 7.3 − 2(6 .5) + 5 .1 12 = −0.6 m / s 2 b) v= a= − x (ti + 2 ) + 4 x (ti +1 ) − 3x (ti ) − 8 + 4(7.3) − 3( 6.5) = = 0.85 m / s 2h 2 − x(ti + 3 ) + 4 x(t i + 2 ) − 5 x(t i +1) + 2 x (ti ) h 2 = − 8.4 + 4(8) − 5(7.3) + 2( 6.3) 2 1 = −0.3 m / s 2 c) v= a= 3 x (ti ) − 4 x (ti −1) + x( ti − 2 ) 3(6 .5) − 4(5.1) + 3.4 = = 1 .25 m / s 2h 2 2 x(t i ) − 5 x(t i −1 ) + 4 x (ti − 2 ) − x (ti − 3 ) h 2 = 2(6.5) − 5(5.1) + 4(3 .4) − 1.8 2 1 = − 0. 7 m / s 23.22 dθ θ (ti +1 ) − θ (ti −1) 0.67 − 0.70 θ = = = = −0.015 rad/s dt 2h 2 r = dr r (ti +1 ) − r (ti −1 ) 6030 − 5560 = = = 235 ft/s dt 2h 2 2 d 2θ θ (ti +1 ) − 2θ (ti ) + θ (t i −1 ) 0 .67 − 2( 0.68) + 0.70 θ = 2 = = = 0 .01 rad/s 2 2 dt h (1) r = d 2r dt 2 = r (ti +1 ) − 2 r (ti ) + r (t i −1 ) h 2 = 6030 − 2 (5800) + 5560 (1) 2 2 = −10 ft/s    v = 235 er − 87 eθ    a = −8.695 er + 50.95 eθ 23.23 Use the same program as was developed in the solution of Prob. 23.11 2 23.24