数值计算

Posted on | In | |
Words count in article: 658 | Reading time ≈ 3

Representing Polynomials in MATLAB

  • Polynomials were represented as row vectors
  • For example: f(x) = x^3 - 2*x - 5
  • to enter this polynomial into MATLAB, use p = [1 0 -1 -5];

values of polynomials : polyval()

  • plot the polynomial

f(x) = 9x^3 - 5x^2 + 3x + 7 (for -2 <= x <= 5)

1
2
3
4
5
6
a = [9, -5, 3, 7];
x = -2:0.01:5;
f = polyval(a, x);
plot(x, f, 'LineWidth', 2);
xlabel('x'); ylabel('f(x)');
set(gca, 'FontSize', 14);

Polynomial Differentiation微分: polyder()

1
2
p = [5 0 -2 0 1];
polyder(p)

算某处微分值：

1
2
p = [5 0  -2 0 1];
polyval(polyder(p), 7)

Polynomial Integration积分 ： polyint()

需要多给matlab一个常熟项，积分后会出现一个常数项

1
2
p = [5 0 -2 0 1];
polyint(p, 3)

在某处的值

1
polyval(polyint(p, 3), 7)

Numerical数值 Differentiation微分

Differences: diff()

  • diff() calculates the differences between adjacent相邻 elements of a vector
1
2
x = [1 2 5 2 1];
diff(x)

obtain the slope坡度 of a line between 2 points (1, 5) and (2, 7)

1
2
x = [1 2]; y = [5 7];
slope = diff(y)./diff(x)

Numerical Differentiation Using diff()

  • Given f(x) = sin(x), find f`(x0) at x0 = pi/2 using h = 0.1
1
2
3
4
x0 = pi/2; h = 0.1;
x = [x0 x0+h];
y = [sin(x0) sin(x0+h)];
m = diff(y)./diff(x)

How to Find the f` over an interval [0, 2Π]?

1
2
3
h = 0.05; x = 0:h:2*pi;
y = sin(x); m = diff(y)./diff(x);
plot(m);

Second and Third Derivatives 二次/三次微分

1
2
3
4
5
6
7
8
9
x = -2:0.005:2; y = x.^3;
m = diff(y)./diff(x);
m2 = diff(m)./diff(x(1:end-1));

plot(x, y, x(1:end-1), m, x(1:end-2), m2);
xlabel('x', 'FontSize', 18);
ylabel('y', 'FontSize', 18);
legend('f(x) = x^3', 'f''(x)', 'f''''(x)');
set(gca, 'FontSize', 18);

1

Numerical Integration

用矩形或梯形来估算面积大小来作为积分值

Midpoint矩形估算 Rule Using sum()

底乘高，底是h ， 高是中等分点的函数值

1
2
3
4
h = 0.05; x = 0:h:2;
midpoint = (x(1:end-1)+x(2:end))./2;
y = 4*midpoint.^3;
s = sum(h*y);

Trapezoid Rule 梯形预测

上底加下底乘高除二

上底是x1的函数值 下底是x0的函数值 高是h

1
2
h = 0.05; x = 0:h:2; y = 4*x.^3;
s = h*trapz(y)  % trapz() 就是计算上底加下底除二

或者：

1
2
3
h = 0.05; x = 0:h:2; y = 4*x.^3;
trapezoid = (y(1:end-1)+y(2:end))/2;
s = h*sum(trapezoid)

Second-order Rule: 1/3 Simpson’s

计算面积的公式为：

h/3*(f0 + 4f1 + f2)

1

Review of Function Handles (@)

  • A handle is a pointer to a function
  • can be used to pass functions to other functions
1
2
3
4
5
function [y] = xy_plot(input, x)
% xy_plot receives the handle of a function and plots that function of x
y = input(x); plot(x, y 'r--');
xlabel('x'); ylabel('function(x)');
end

Try: xy_plot(@sin, 0:0.01:2*pi);

不能将函数传给函数，只能传一个函数的handle给另一个函数

Numerical Integration积分: integral()

1
2
y = @(x) 1./(x.^3-2*x-5);
integral(y,0,2)

Double and Triple Integrals

1

Donate? comment?