課件文稿成果lecture

1、and its application inEngineeringAssoc. Prof. Kirin Shi施圣賢上海交通大學機械與動力課程的作業(yè)上傳至如下ftp中的upload文件夾ftp地址:（在ftp:/管理器的地址欄中輸入，不是瀏覽器）用戶名: 0120209060: public課程助教:華晨Review on Lecture 3 Symbolic match Differentiation Integration Representation of polynomials Solving polynomial equationSymbolic mathDeclaring symb

2、olic variablessyms symbolic variables options syms x a real>> f = (x-1)*(x-a)*(2*x2+5)+ 20syms without a list of variables, you would get the list of all the symbolic variables in your workspace>> syms'a''f''x'Symbolic mathPerform symbolic calculaton>> f_exp

3、and = expand(f)f_expand =2*x4+5*x2-2*x3*a-5*x*a-2*x3-5*x+2*x2*a+5*a+20>> a=2;>> f_new = eval(f)>> f_expand_new = eval(f_expand) f_expand_new =2*x4+9*x2-6*x3-15*x+30>> expand(sin(2*x)Symbolic mathPerform symbolic calculationuse the function simple to find the simplest (shortes

4、t)form of a symbolic expression>> R = simple(sin(x)2 + cos(x)2) R =1>> R = simple(x*(x+1)*(x+2)>> R = simple(x3+3*x2+3*x+1)>> R = simple(2*sin(x)*cos(x)+2)Differentiationfunctiondiff(y,x,n)calculate the 1st and 2ndderivative of a functiony(x) = sin(a*x) + 4x5 with respect to

5、the variable x>> syms a x>> y_first_derivative = diff(y,x)>> y_second_derivative = diff(y,x,2)>> y_second_derivative_a = diff(y,a,2)Integrationfunctionint(f,x)int(f,x,a,b)bA = ò f ( x)dxaA = ò f ( x)dx¥>> syms a x-2|x|A =òedx-¥>> f = cos(a

6、*x)*a+20*x4>> syms a b x>> f = x2 + 3*x +1;A = int(exp(-2*abs(x),-inf,inf)>> B = int(f,x,a,b)Integrationtrapz, quad and dblquadz = trapz (x,y) computes the integral of y with respect to x using the trapezoidal rule3 sin 2 (x)A = òdxx1X = linspace (1,3,11); Y = (sin(x).2)./x;A

7、= trapz (x, y)A = 0.7938Effect of more intervalsx=linspace(1,3,21);y=(sin(x).2)./x;x=linspace(1,3,41);y=(sin(x).2)./x;A = trapz(x,y)A = 0.7946A = trapz(x,y)A = 0.7948>> syms x>> Y = (sin(x).2)./x;>> int(Y,x,1,3)>> cosint(2)/2 - cosint(6)/2 + log(3)/2ans =0.7948Function quadNu

8、merically evaluating the integral of a function using adaptive Simpson quadrature, to within an error of 1.e-6.Q = quad (FUN, lower_limit, upper_limit)A = quad ('sin(x).2./x', 1, 3)A = 0.7948B = quad ('exp(-x.2)', 0, 1)B = 0.7468sin2 (x)3A = òdx;x11B = òe-x2 dx0Function dbl

9、quadnumerically evaluate a double integral using adaptive Simpson quadrature.15A = ò ò x2 y3 dx dy-1 -1Q = dblquad(FUN, x1, x2, y1, y2)A = dblquad ('x.2.*y.3', -1, 1, -1, 5) A =104.0000quadgk: adaptive Gauss-Kronrod quadraturequad2d: Numerically evaluate integral, adaptive Lobatto

10、quadraturePOLYNOMIALSMany physical and engineering problems lead to polynomialfunctions.A polynomial function of degree n has teral formy = an xn + an-1 xn-1 + . + a1 x + a0A solution or root of the equation is a number xi, real or complex, such thatan xn + an-1 xn-1 + . + a1 xi + a0 = 0iiRepresenta

11、tion of polynomialsRepresented by the vector of its coefficients, in descending powers of x.The length of the coefficient vector is always equal to theorder of the polynomial plus one.For example,x2 +x +x x-xFor example,>> coef_poly_1 = 1coef_poly_1 =21+ 2x + 1x2121>> coef_poly_2 = 2coef

12、_poly_2 =0-15-12- 15x -122x320-15-12>> coef_poly_3 = 2coef_poly_3 =0-150- 15x2x320-150Solving polynomial equationsUse the function roots for finding the roots of a polynomial equation.The argument of this function is the vector of the polynomial coefficients.+ 2x +1x2>> coef_poly_1 = 121

13、>> roots_eq_1 = roots(coef_poly_1)roots_eq_1 =-1-1Function fzeroIn general, we would use the function fzero forsolving equations.The commandfzero ('fun', X0) tries to find azero (a root) of the function fun near X0.equation = '2*x.4-55*x.3 + 506.02*x.2 - 1755.89*x + 1670.79'R1

14、= fzero (equation, -1)R1 = 1.5000>> R1 = fzero (equation, 3)R1 =1.5000>> R1 = fzero (equation, 4)R1 =5.5000>> R1 = fzero (equation, 7)R1 =5.5000>> R1 = fzero (equation, 8)R1 =8.3000>> R1 = fzero (equation, 12)R1 =12.2000>> R1 = fzero (equation, 16)R1 =12.2000y = 2

15、x4 - 55x3 + 506.02x2 - 1755.89x + 1670.7925002000150010005000-50002468101214xyMultiplication and division of polynomials- 2) = ?(2>> coef_poly_1 = 2 3 1;>> coef_poly_2 = 5 -2;>> coef_product = conv(coef_poly_1,coef_poly_2);>> coef_product = coef_product(:) coef_product =10.00

16、11.00-1.00-2.00- 210Multiplication and division of polynomials+ 2 =+ 5x - 28?x2>> coef_poly_num = 8 12 -3 2;>> coef_poly_den = 1 5 -2;>>coef_quotient,coef_remainder=deconv(coef_poly_num,coef_poly_den)The first one is the coefficients of the quotient polynomialThe second one is the

17、coefficients of the remainder polynomialcoef_quotient = 8.00coef_remainder =0-28.000153.00-54.008x3 + 12+ 5- 54)2Non-polynomial functionsin( x) = ex -5sin( x) - ex +5 = 0root1 = fzero ('sin(x)-exp(x)+5', 16)root1 = 1.7878function f1200-20-40-60-80-100-120-140-160-5-4-3-2-10x12345yNon-polynom

18、ial function: Define equation by m-filesin( x) = ex -5sin( x) - ex +5 = 0%function M-file to define f1 function y = f1(x)y = sin(x) - exp(x) + 5;>> root1 = fzero('f1',16)root1 = 1.7878function f1200-20-40-60-80-100-120-140-160-5-4-3-2-10x12345yDefining a polynomial with known rootsThe

19、function poly is the inverse of the function roots.For example, determine the coefficients of a polynomial whose roots are equal to1+i1+i 2+3i4>> equation_roots>> equation_coefs>> equation_coefs equation_coefs = 1.00002.0000 - 5.0000i-16.0000 - 16.0000i-28.0000 + 22.0000i16.0000 +

20、24.0000i=-1+i,1+i,2+3i,-4;poly(equation_roots);equation_coefs(:)x4 + (2 - 5i) x3 - (16 +16i) x2 + (-28 + 22i) x + (16 + 24i)Finding the values of a polynomialOrdinary procedureFor example, the values of the functiony = x5 - 2+ 42for x = -2:0.1:2>> x = -2:0.1:2;>> y = x.5 -2*x.4 + 2*x.3 +

21、3*x.2 + 4;function polyvaly = x5 - 2+ 42>> coef_poly_3 = 1-22304;>> x = -2:0.1:2;>> y = polyval(coef_poly_3,x);Fitting a polynomial to dataxyxy-2.00-8.430.604.48-1.80-6.380.801.50-1.60-1.891.001.11-1.400.421.201.51-1.200.651.400.08-1.004.191.60-0.32-0.805.001.80-1.08Variation of y

22、with x from laboratory data86420-2-4-6-8-10-2-10123xVariation of y with x from experimental datayFitting a polynomial of order n to the data points (xi, yi)>> coef_poly = polyfit(x,y,n)n is the order of the polynomial.coef_poly is the vector of the coefficients of a polynomial of order n that

23、best fits the data (xi, yi) in a least square sense.We can then use the vector coef_poly to evaluate thepolynomial for the given x values>> x = -2, -1.8, .;>> y = -8.43, -6.38, .;>> coef_poly = polyfit(x,y,3);>> y_hat = polyval(coef_poly, x);>> figure, plot(x, y_hat, &#

24、39;-', x, y, 'o')Third order polynomial fitted to experimental data86420-2Fitted polynomialObserved data-4-6-8-10-2-10123xThird order polynomial fitted to experimental data yCURVE FITTING by LSETechnique used for fitting an algebraic function to data.Usually done by the Least Square Erro

25、r (LSE) method.xyxyxy0.0030.997.0018.9514.0036.030.5022.677.5022.8314.5035.251.0018.968.0023.5015.0032.991.5018.178.5022.5615.5024.942.0015.629.0020.0216.0018.852.5016.149.5016.7316.5010.463.0015.7810.0010.8417.007.793.5011.7910.507.3917.506.224.0012.0111.005.7818.0010.744.509.0811.509.7318.5019.315

26、.008.2812.0011.7819.0026.795.5010.7612.5018.1219.5036.836.0012.2813.0025.4120.0045.176.5015.0813.5031.28Graphical representation of the experimental dataTheoretical or empirical relationshipay =+ b x+ c xsin( x) + dx +1ayi =+ b xi+ c xi sin( xi ) + dxi + 1i = 1, 2,., 41Derivation of parametersOnly 4

27、 unknown parameters (a, b, c and d), but41 equations.Equations are said to be overdetermined.In other words, no set of a, b, c and d values would satisfy all the 41 equations.Consequently, for any choice of the parameters, the 41 equations can be written in the following formæöa- çx+

28、c xi sin( xi ) + d ÷ = ei+ 1 + b xiyièøii = 1, 2,., 41LSE Methodei is the experimental error and/or inaccuracies inthe theoretical m.Determine the parameters a, b, c and d to minimise the error terms ei in some sense. not good to minimise the sum of the errors as the positive and nega

29、tive errors can cancel each other out.Minimise the sum of the squares of the errors, because the squares are always positive and they do not cancel each other out.a+ bx+ cx sin(x ) + d = yx + 1iiiiii = 1, 2,., 41é1ù1xx sin(x )ê xúúú+ 1111êê1é y1ù

30、0;úúúúúúúú1xx2 sin(x2 ).1úú.ú.úê yê x+ 1. 12êêêêêêêê22éaù.êbúêúúê úêcú=ê yixx sin(x )1úú.ú.úêêê

31、 xêdú+ 1.iii.ë ûêêêiê.êyúêë41 úûêúú1úûê1ê xxx41 sin(x41 )+ 141ë41= y A * parametersDefine matrix A and the vertical vector y.Determine parameters a, b, c and d fromparameters = Ay.Always use

32、 the backlash operator to solve a set of linear equations.If the number of unknowns and equations are equal to each other, there would be a unique set of parameter values which would satisfy allthe equations. = yA * parametersIf the number of equations are more than the number of parameters, then th

33、e equations will be solved using the Least Square Error technique.Commandsx = 0:0.5:20;y =30.99 22.67 18.96 18.17 15.62 16.14 15.7811.7922.8312.01 9.088.28 10.76 12.28 15.08 18.9523.50 22.56 20.02 16.73 10.84 7.39 5.789.73 11.78 18.12 25.41 31.28 36.03 35.25 32.9924.94 18.85 10.46 7.79 6.22 10.74 19

34、.31 26.7936.83 45.17;x = x(:);y = y(:);N = length(x);= yA * parameters>> A = 1./(x+1)xx.*sin(x)ones(N,1);>> C = AyPlot a smooth curve through the data pointsxn = linspace (0, 20, 200);y_p = C(1)./(xn+1)+C(2)*xn + . C(3)*xn.*sin(xn)+C(4);figure, plot (x, y, 'o', xn, y_p, '-

35、9;), grid on, axis squarexlabel('x','Fontsize',14),ylabel('y','Fontsize',14) title('Fitting curve ', 'Fontsize',14)legend('Measured data',' LSE')50INTERPOLATIONHave a set of x and y values.Pass a smooth curve through these pointsor find

36、 the value of y for some intermediate xvalues.This operation is called interpolation.Not same as curve fitting, as in curve fitting, the fitted curve does not, in general, pass through the data points.function interp1The function interp1 can be used for one-dimensional interpolation.Teral form of th

37、is function isy_new = interp1 (x, y, x_new, 'method')nearest, linear, 'spline, cubiccommands>> x_new = linspace(0,20,200);>> y_new = interp1(x,y, x_new,'spline');>> figure, plot (x, y, 'o', x_new, y_new, '-')>> grid on, axis square54functio

38、n interp2Two-dimensional interpolation.Teral form of this function isy_new = interp2 (x, y, f(x,y), x_new, y_new'method')nearest, linear, 'spline, cubiccommandsX,Y = meshgrid(-3:0.25:3);Z = peaks(X,Y);XI,YI = meshgrid(-3:0.125:3);ZI = interp2(X,Y,Z,XI,YI);figure, mesh(X,Y,Z), hold on, me

39、sh(XI,YI,ZI+15)hold off, axis(-3 3 -3 3 -5 20)FINDING THE MINIMUMOF A FUNCTIONy = 2x4 - 55x3 + 506.02x2 - 1755.89x + 1670.7925002000150010005000-50002468101214xyThe command fminsearchFminsearch ('fun', X0) tries to find a local minimum of the function fun near X0.>> equation = '2*x

40、.4-55*x.3+506.02*x.2-1755.89*x.79'>> x_min_1 = fminsearch (equation,12)x_min_1 = 10.7745>> x_min_2 = fminsearch (equation, 5)x_min_2 = 2.9537Minimum of a function defined by M-filefunction y = f_poly_1(x)y = 2*x.4 - 55*x.3 + 506.02*x.2 -1755.89*x + 1670.79;>> x_min_1 = fminsear

41、ch('f_poly_1',12)x_min_1 = 10.7745>> x_min_2 = fminsearch('f_poly_1',5)x_min_2 = 2.9537The value of the function at its minimum>> y_min_1 = f_poly_1 (x_min_1)y_min_1 =-345.1092>> y_min_2 = f_poly_1(x_min_2)y_min_2 =-365.9617Quick Excise1) Write a function M-file to

42、calculate y with given input x accordingto𝑦 = 3𝑒𝑥 + sin 4𝑥 𝑥 𝑥2) Then fits a polynomial of order 5, 10, to the (x, y) data.3) Plot y, and the fitted polynomial against the given x values.4) Interpolate y (y_new) according to new x values (x_new) by usi

43、ng spline interpolation method5) Plot x, yand x_new, y_newFunction inputs: x=-2:0.2:3, x_new= -1:0.05:2Functio outputs: figure 1,figure 2,figure 3,plot x,y and x5,y5; plot x, y and x10,y10plot x, y and x_new, y_newHome work基礎及其應用Section 5.7 P154Section 6.9 P183Section 7.6 P212STATISTICALANALYSIS AND

44、 SIMULATIONRandom Samples from aUniform Distribution between zero and one>> xx = rand(3,5)xx = 0.67210.83810.01960.68130.37950.83180.50280.70950.42890.30460.18970.19340.68220.30280.5417Random Samples from a Uniform Distribution between a and b>> a = 10;>> b = 85;>> xx = rand(

45、3,5);>> yy = (b-a)*xx + ayy =60.410372.858911.473061.095838.461172.384747.711063.210442.166932.846324.224024.507361.166732.707350.6255Random Integers between -20 & 150First simulate random numbers between -20 and 150>> a = -20;>> b = 150;>> xx = rand(3,5);>> yy = (b

46、-a)*xx+ayy =72.6372121.2067109.0900126.7859123.1127-12.1738100.3727102.597723.815615.742888.2388-4.107791.0830136.046477.3836Then use the function round to convert the realnumbers into integers.>> zz = round(yy)zz =73121109127123-12100103241688-49113677Standardised Gaussian (Normal) random var

47、iableStandardised Gaussian random variable0.4E(X) = 0;(X) = 1;0.350-4-20x24probability density functionRandom Samples from a Standardised Gaussian Distribution>> xx = randn(3,5)xx = 2.02110.5018-1.99830.27230.33680.1378-1.6106-1.0075-0.5144-2.08891.04610.2153-0.1624-0.1758

48、1.1022Gaussian random variableGaussian random variable, meany=175, stdy = 200.02mean_y = 175;sigma_y = 20;0.0150.010.005050100150200250300yprobability density functionRandom Samples from a Gaussian Distribution>> mean_y = 175;>> sigma_y = 20;First simulate random numbers from a standardi

49、sed Gaussian random variable.>> x = randn(3,5)x =0.89790.50580.3299-0.09010.64811.7347-0.0858-0.30831.4596-1.4038-0.70920.63300.2215-1.33260.7305Random Samples from a Gaussian DistributionThen use the following relationship to change the standard deviation and the mean of the random variable.&

50、gt;> y = sigma_y*x + mean_yy =192.9581185.1164181.5986173.1973187.9621209.6938173.2833168.8337204.1930146.9244160.8168187.6603179.4303148.3489189.6097Sample mean and standard deviation>> mean_y = 175;>> sigma_y = 20;>> x = randn(1,1000);mean_y =174.7924>> y = sigma_y*x + m

51、ean_y;std_y =19.6464>> mean_y = mean(y)>> std_y = std(y)Uniform distribution>> a=1000;>> b=10000;>> xx = rand(1,2500);>> yy = (b-a)*xx + a;>> hist(yy,20)>> xlabel('bins')>> ylabel('number of elements in each bin')>> title('Histogram