[Matlab-Based Book]MATLAB程式設計與應用(第三版)

本书将 MATLAB 作为一种技术编程语言进行介绍，并且向学生展示如何编写效率高、组织严谨的程序，同时强调问题解决技巧。书中重点介绍良好的编程实践、MATLAB 编程中的常见错误，以及大量使用 MATLAB 科学计算软件解决的练习和示例。

1. %  Script file: calc_roots.m
   
2. %
   
3. %  Purpose:
   
4. %    This program solves for the roots of a quadratic equation
   
5. %    of the form a*x**2 + b*x + c = 0.  It calculates the answers
   
6. %    regardless of the type of roots that the equation possesses.
   
7. %
   
8. %  Record of revisions:
   
9. %      Date       Programmer          Description of change
   
10. %      ====       ==========          =====================
    
11. %    01/12/07    S. J. Chapman        Original code
    
12. %
    
13. % Define variables:
    
14. %   a            -- Coefficient of x^2 term of equation
    
15. %   b            -- Coefficient of x term of equation
    
16. %   c            -- Constant term of equation
    
17. %   discriminant -- Discriminant of the equation
    
18. %   imag_part    -- Imag part of equation (for complex roots)
    
19. %   real_part    -- Real part of equation (for complex roots)
    
20. %   x1           -- First solution of equation (for real roots)
    
21. %   x2           -- Second solution of equation (for real roots)
    
22. 
    
23. % Prompt the user for the coefficients of the equation
    
24. disp ('This program solves for the roots of a quadratic ');
    
25. disp ('equation of the form A*X^2 + B*X + C = 0. ');
    
26. a = input ('Enter the coefficient A: ');
    
27. b = input ('Enter the coefficient B: ');
    
28. c = input ('Enter the coefficient C: ');
    
29. 
    
30. % Calculate discriminant
    
31. discriminant = b^2 - 4 * a * c;
    
32. 
    
33. % Solve for the roots, depending on the value of the discriminant
    
34. if discriminant > 0 % there are two real roots, so...
    
35. 
    
36.    x1 = ( -b + sqrt(discriminant) ) / ( 2 * a );
    
37.    x2 = ( -b - sqrt(discriminant) ) / ( 2 * a );
    
38.    disp ('This equation has two real roots:');
    
39.    fprintf ('x1 = %f\n', x1);
    
40.    fprintf ('x2 = %f\n', x2);
    
41. 
    
42. elseif discriminant == 0  % there is one repeated root, so...
    
43. 
    
44.    x1 = ( -b ) / ( 2 * a );
    
45.    disp ('This equation has two identical real roots:');
    
46.    fprintf ('x1 = x2 = %f\n', x1);
    
47. 
    
48. else % there are complex roots, so ...
    
49. 
    
50.    real_part = ( -b ) / ( 2 * a );
    
51.    imag_part = sqrt ( abs ( discriminant ) ) / ( 2 * a );
    
52.    disp ('This equation has complex roots:');
    
53.    fprintf('x1 = %f +i %f\n', real_part, imag_part );
    
54.    fprintf('x1 = %f -i %f\n', real_part, imag_part );
    
55. 
    
56. end

复制代码

1. %  Script file: funxy.m
   
2. %
   
3. %  Purpose:
   
4. %    This program solves the function f(x,y) for a
   
5. %    user-specified x and y, where f(x,y) is defined as:
   
6. %                 _
   
7. %                |
   
8. %                | x + y             x >= 0 and y >= 0
   
9. %                | x + y^2           x >= 0 and y < 0
   
10. %       f(x,y) = | x^2 + y           x < 0  and y >= 0
    
11. %                | x^2 + y^2         x < 0  and y < 0
    
12. %                |_
    
13. %
    
14. %  Record of revisions:
    
15. %      Date       Programmer          Description of change
    
16. %      ====       ==========          =====================
    
17. %    01/13/07    S. J. Chapman        Original code
    
18. %
    
19. % Define variables:
    
20. %   x     -- First independent variable
    
21. %   y     -- Second independent variable
    
22. %   fun   -- Resulting function
    
23. 
    
24. % Prompt the user for the values x and y
    
25. x = input ('Enter the x coefficient: ');
    
26. y = input ('Enter the y coefficient: ');
    
27. 
    
28. % Calculate the function f(x,y) based upon
    
29. % the signs of x and y.
    
30. if x >= 0 && y >= 0
    
31.    fun = x + y;
    
32. elseif x >= 0 && y < 0
    
33.    fun = x + y^2;
    
34. elseif x < 0 && y >= 0
    
35.    fun = x^2 + y;
    
36. else
    
37.    fun = x^2 + y^2;
    
38. end
    
39. 
    
40. % Write the value of the function.
    
41. disp (['The value of the function is ' num2str(fun)]);

复制代码

1. %  Script file: ideal_gas.m
   
2. %
   
3. %  Purpose:
   
4. %    This program plots the pressure versus volume of an
   
5. %    ideal gas.
   
6. %
   
7. %  Record of revisions:
   
8. %      Date       Programmer          Description of change
   
9. %      ====       ==========          =====================
   
10. %    01/16/07    S. J. Chapman        Original code
    
11. %
    
12. % Define variables:
    
13. %   n         -- Number of atoms (mol)
    
14. %   P         -- Pressure (kPa)
    
15. %   R         -- Ideal gas constant (L kPa/mol K)
    
16. %   T         -- Temperature (K)
    
17. %   V         -- volume (L)
    
18. 
    
19. % Initialize nRT
    
20. n = 1;                   % Moles of atoms
    
21. R = 8.314;               % Ideal gas constant
    
22. T = 273;                 % Temperature (K)
    
23. 
    
24. % Create array of input pressures.  Note that this
    
25. % array must be quite dense to catch the major
    
26. % changes in volume at low pressures.
    
27. P = 1:0.1:1000;
    
28. 
    
29. % Calculate volumes
    
30. V = (n * R * T) ./ P;
    
31. 
    
32. % Create first plot
    
33. figure(1);
    
34. loglog( P, V, 'r-', 'LineWidth', 2 );
    
35. title('\bfVolume vs Pressure in an Ideal Gas');
    
36. xlabel('\bfPressure (kPa)');
    
37. ylabel('\bfVolume (L)');
    
38. grid on;
    
39. hold on;
    
40. 
    
41. % Now increase temperature
    
42. T = 373;                 % Temperature (K)
    
43. 
    
44. % Calculate volumes
    
45. V = (n * R * T) ./ P;
    
46. 
    
47. % Add second line to plot
    
48. figure(1);
    
49. loglog( P, V, 'b--', 'LineWidth', 2 );
    
50. hold off;
    
51. 
    
52. % Add legend
    
53. legend('T = 273 K','T = 373 K');

复制代码

1. %  Script file: ball.m
   
2. %
   
3. %  Purpose:
   
4. %    This program calculates the distance traveled by a ball
   
5. %    thrown at a specified angle "theta" and a specified
   
6. %    velocity "vo" from a point on the surface of the Earth,
   
7. %    ignoring air friction and the Earth's curvature.  It
   
8. %    calculates the angle yielding maximum range, and also
   
9. %    plots selected trajectories.
   
10. %
    
11. %  Record of revisions:
    
12. %      Date       Programmer          Description of change
    
13. %      ====       ==========          =====================
    
14. %    01/30/07    S. J. Chapman        Original code
    
15. %
    
16. % Define variables:
    
17. %   conv         -- Degrees to radians conv factor
    
18. %   gravity      -- Accel. due to gravity (m/s^2)
    
19. %   ii, jj       -- Loop index
    
20. %   index        -- Location of maximum range in array
    
21. %   maxangle     -- Angle that gives maximum range (deg)
    
22. %   maxrange     -- Maximum range (m)
    
23. %   range        -- Range for a particular angle (m)
    
24. %   time         -- Time (s)
    
25. %   theta        -- Initial angle (deg)
    
26. %   traj_time    -- Total trajectory time (s)
    
27. %   vo           -- Initial velocity (m/s)
    
28. %   vxo          -- X-component of initial velocity (m/s)
    
29. %   vyo          -- Y-component of initial velocity (m/s)
    
30. %   x            -- X-position of ball (m)
    
31. %   y            -- Y-position of ball (m)
    
32. 
    
33. %  Constants
    
34. conv = pi / 180;      % Degrees-to-radians conversion factor
    
35. g = -9.81;            % Accel. due to gravity
    
36. vo = 20;              % Initial velocity
    
37. 
    
38. %Create an array to hold ranges
    
39. range = zeros(1,91);
    
40. 
    
41. % Calculate maximum ranges
    
42. for ii = 1:91
    
43.    theta = ii - 1;
    
44.    vxo = vo * cos(theta*conv);
    
45.    vyo = vo * sin(theta*conv);
    
46.    max_time = -2 * vyo / g;
    
47.    range(ii) = vxo * max_time;
    
48. end
    
49. 
    
50. % Write out table of ranges
    
51. fprintf ('Range versus angle theta:\n');
    
52. for ii = 1:91
    
53.    theta = ii - 1;
    
54.    fprintf('  %2d    %8.4f\n',theta, range(ii));
    
55. end
    
56. 
    
57. % Calculate the maximum range and angle
    
58. [maxrange index] = max(range);
    
59. maxangle = index - 1;
    
60. fprintf ('\nMax range is %8.4f at %2d degrees.\n',...
    
61.          maxrange, maxangle);
    
62. 
    
63. % Now plot the trajectories
    
64. for ii = 5:10:85
    
65. 
    
66.    % Get velocities and max time for this angle
    
67.    theta = ii;
    
68.    vxo = vo * cos(theta*conv);
    
69.    vyo = vo * sin(theta*conv);
    
70.    max_time = -2 * vyo / g;
    
71. 
    
72.    % Calculate the (x,y) positions
    
73.    x = zeros(1,21);
    
74.    y = zeros(1,21);
    
75.    for jj = 1:21
    
76.       time = (jj-1) * max_time/20;
    
77.       x(jj) = vxo * time;
    
78.       y(jj) = vyo * time + 0.5 * g * time^2;
    
79.    end
    
80.    plot(x,y,'b');
    
81.    if ii == 5
    
82.       hold on;
    
83.    end
    
84. end
    
85. 
    
86. % Add titles and axis lables
    
87. title ('\bfTrajectory of Ball vs Initial Angle \theta');
    
88. xlabel ('\bf\itx \rm\bf(meters)');
    
89. ylabel ('\bf\ity \rm\bf(meters)');
    
90. axis ([0 45 0 25]);
    
91. grid on;
    
92. 
    
93. % Now plot the max range trajectory
    
94. vxo = vo * cos(maxangle*conv);
    
95. vyo = vo * sin(maxangle*conv);
    
96. max_time = -2 * vyo / g;
    
97. 
    
98. % Calculate the (x,y) positions
    
99. x = zeros(1,21);
    
100. y = zeros(1,21);
     
101. for jj = 1:21
     
102.    time = (jj-1) * max_time/20;
     
103.    x(jj) = vxo * time;
     
104.    y(jj) = vyo * time + 0.5 * g * time^2;
     
105. end
     
106. plot(x,y,'r','LineWidth',3.0);
     
107. hold off

复制代码

1. %  Script file: doy.m
   
2. %
   
3. %  Purpose:
   
4. %    This program calculates the day of year corresponding
   
5. %    to a specified date.  It illustrates the use switch
   
6. %    and for constructs.
   
7. %
   
8. %  Record of revisions:
   
9. %      Date       Programmer          Description of change
   
10. %      ====       ==========          =====================
    
11. %    01/27/07    S. J. Chapman        Original code
    
12. %
    
13. % Define variables:
    
14. %   day          -- Day (dd)
    
15. %   day_of_year  -- Day of year
    
16. %   ii           -- Loop index
    
17. %   leap_day     -- Extra day for leap year
    
18. %   month        -- Month (mm)
    
19. %   year         -- Year (yyyy)
    
20. 
    
21. % Get day, month, and year to convert
    
22. disp('This program calculates the day of year given the ');
    
23. disp('specified date.');
    
24. month = input('Enter specified month (1-12): ');
    
25. day   = input('Enter specified day(1-31):    ');
    
26. year  = input('Enter specified year(yyyy):   ');
    
27. 
    
28. % Check for leap year, and add extra day if necessary
    
29. if mod(year,400) == 0
    
30.    leap_day = 1;          % Years divisible by 400 are leap years
    
31. elseif mod(year,100) == 0
    
32.    leap_day = 0;          % Other centuries are not leap years
    
33. elseif mod(year,4) == 0
    
34.    leap_day = 1;          % Otherwise every 4th year is a leap year
    
35. else
    
36.    leap_day = 0;          % Other years are not leap years
    
37. end
    
38. 
    
39. % Calculate day of year by adding current day to the
    
40. % days in previous months.
    
41. day_of_year = day;
    
42. for ii = 1:month-1
    
43. 
    
44.    % Add days in months from January to last month
    
45.    switch (ii)
    
46.    case {1,3,5,7,8,10,12},
    
47.       day_of_year = day_of_year + 31;
    
48.    case {4,6,9,11},
    
49.       day_of_year = day_of_year + 30;
    
50.    case 2,
    
51.       day_of_year = day_of_year + 28 + leap_day;
    
52.    end
    
53. 
    
54. end
    
55. 
    
56. % Tell user
    
57. fprintf('The date %2d/%2d/%4d is day of year %d.\n', ...
    
58.          month, day, year, day_of_year);

复制代码

1. %  Script file: logical1.m
   
2. %
   
3. %  Purpose:
   
4. %    This program calculates the time required to
   
5. %    calculate the square roots of all elements in
   
6. %    array a whose value exceeds 5000.  This is done
   
7. %    in two different ways:
   
8. %    1.  Using a for loop and if construct.
   
9. %    2.  Using a logical array.
   
10. %
    
11. %  Record of revisions:
    
12. %      Date       Programmer          Description of change
    
13. %      ====       ==========          =====================
    
14. %    01/27/07    S. J. Chapman        Original code
    
15. %
    
16. % Define variables:
    
17. %   a            -- Array of input values
    
18. %   b            -- Logical array to serve as a mask
    
19. %   ii, jj       -- Loop index
    
20. %   average1     -- Average time for calculation 1
    
21. %   average2     -- Average time for calculation 2
    
22. %   maxcount     -- Number of times to loop calculation
    
23. %   month        -- Month (mm)
    
24. %   year         -- Year (yyyy)
    
25. 
    
26. 
    
27. % Perform calculation using loops and branches.
    
28. maxcount = 10;              % Number of repetitions
    
29. tic;                        % Start timer
    
30. for jj = 1:maxcount        
    
31.    a = 1:10000;             % Declare array a
    
32.    for ii = 1:10000   
    
33.       if a(ii) > 5000
    
34.          a(ii) = sqrt(a(ii));
    
35.       end
    
36.    end
    
37. end
    
38. average1 = (toc)/maxcount;  % Calculate average time
    
39. 
    
40. 
    
41. % Perform calculation using logical arrays.
    
42. maxcount = 100;             % Number of repetitions
    
43. tic;                        % Start timer
    
44. for jj = 1:maxcount        
    
45.    a = 1:10000;             % Declare array a
    
46.    b = a > 5000;            % Create mask
    
47.    a(b) = sqrt(a(b));       % Take square root
    
48. end
    
49. average2 = (toc)/maxcount;  % Calculate average time
    
50. 
    
51. % Display results
    
52. fprintf('Loop / if approach =     %8.4f\n', average1);
    
53. fprintf('Logical array approach = %8.4f\n', average2);
    
54.          

复制代码