[转帖]MATLAB多项式两个不错的例子

aspen

来自：http://www.rit.edu/~pnveme/pigf/Polynomials/poly_index.html

例一：问题描述




定义：



For w = 60,000 N/m; L = 6 m; E = 200 GPa; I = 15E-06 kg - m2

1. % in MATLAB
   
2. % polynomial described as a vector of coefficients
   
3. 
   
4. p = [0.000833 -0.01 0 0.18 0]
   
5. 
   
6. % roots of the polynomial - points where ploynomial is zero
   
7. % we know that it is zero at 0 and 6
   
8. % (but it has 4 roots)
   
9. r = roots(p) % roots are in r
   
10. 
    
11. % we can obtain the original polynomial
    
12. % by using the poly function
    
13. p1 = poly(r)
    
14. 
    
15. % returns a normalized polynomial
    
16. % poly and roots are inverse functions
    
17. 
    
18. % symbolic computation - expresses the polynomial
    
19. % from the numerical format - note the rational
    
20. % format for coefficients
    
21. p2 = poly2sym(p)
    
22. 
    
23. % expresses the polynomial numerically
    
24. p3 = sym2poly(p2)
    
25. 
    
26. % change the variable to s
    
27. p4 = poly2sym(p3,'s')
    
28. 
    
29. % in MATLAB
    
30. % polynomial described as a vector of coefficients
    
31. 
    
32. p = [0.000833 -0.01 0 0.18 0];
    
33. 
    
34. % evaluate the polynomial between x = 0 and x = 6
    
35. % and plot the resulting derflaction
    
36. 
    
37. x = 0:0.1:6;
    
38. y = polyval(p,x);
    
39. 
    
40. % since usually the deflection is downwards we plot -y
    
41. plot(x,-y,'r-')
    
42. title('Uniformly Loaded Beam - deflection')
    
43. xlabel('Length along beam - meter')
    
44. ylabel('Deflection of beam - meter')
    
45. text(1,0,'w=60,000, L = 6, E = 200e-09, I = 15 e-06')
    
46. 
    
47. % Find maximum deflection by setting the
    
48. % derivative of polynomial to zero
    
49. q = polyder(p);
    
50. xmax = roots(q);
    
51. for i=1:length(xmax)
    
52. if isreal(xmax(i)) & (xmax(i) > 0) & (xmax(i) < 6)
    
53. ymax=polyval(p,xmax(i));
    
54. fprintf('location of maximum deflection :'),disp(xmax(i))
    
55. fprintf('value of maximum deflection :'),disp(ymax)
    
56. end
    
57. end
    
58. 曲线拟合
    
59. 
    
60. Polynomials are very popular in curve fitting and estimation. We will generate a set of x and y data and will fit the data with different polynomial of a chosen order and compare the approximation
    
61. 
    
62. % in MATLAB
    
63. % create sample data
    
64. x = 0:0.1:3;
    
65. y = sin(2*x);
    
66. % lets fit curves 3,4,5 order
    
67. plot(x,y,'ro','MarkerFaceColor','y')
    
68. hold on
    
69. str={'r--','b-','k:'};
    
70. for i = 3:5
    
71. ii = i -2;
    
72. p = polyfit(x,y,i);
    
73. yy=polyval(p,x);
    
74. plot(x,yy,char(str(ii)));
    
75. legend('Original data','degree 3','degree 4','degree 5')
    
76. end
    
77. hold off
    
78. title('Polynomial Curve fit')
    
79. xlabel('x')
    
80. ylabel('y')
    
81. grid

复制代码

aspen

例二：质量、弹簧、减震系统


The transfer function based on the differential equation can be written in Matrix form as:

The left hand side refers to the Response variables while the right hand side is the Input variables

The above Matrix expression can be written succintly as:

Characteistic Equation:
The characteristic equation is wriiten as:

det A = 0

The roots of this equation indicates if the system is stable or unstable

Using appropriate values for the system parameters the characteristic matrix A can be obtained as



To find the roots of the characteristic equation we have to set the determinant of A to be zero.
The determinant is obtained by the operation (p1 * p3 - p2*p2)
Polynomial multiplication is termed as convolution

    % in MATLAB
    % create the three polynomials
    p1 = [1 6.88 11.8];
    p2 = [-1.6 -14.8];
    p3 = [8.88 58.24 124];

    % you can only add and subtract polynomials
    % of the same order

    p13 = conv(p1,p3);
    p22 = conv(p2,p2);

    % the degree of p13 is 2 more than p22
    % we introduce two leading coefficients of 0]

    p22 = [0 0 p22];
    det = p13 - p22
    % roots of the polynomial
    r = roots(det)

star198311

第一个中的p = [0.000833 -0.01 0 0.18 0]应该是p = [0.000833 0 -0.01 0.18 0]吧？