10

# 數值微分

### 微分

Python:

``````# 函數微分 中央分差
def num_diff(f, x):
h = 1e-4 # 0.0001
return (f(x + h) - f(x - h)) / (2 * h)
``````

Python:

``````# 函數1
def fun_1(x):
return x ** 2 + x

# 微分切線方程式: y = f(a) + d(x − a) = f(a) - d * a + d * x
# 這部分回傳一個function讓我們可以得到函數再把x帶入.
def tangent_fun(f, a):
d = num_diff(f, a)
#print(d)
y = f(a) - d * a
return lambda x: d * x + y

x = np.arange(0.0, 20, 00.1)
y = fun_1(x)

tan_fun = tangent_fun(fun_1, 5)
y2 = tan_fun(x)

plt.plot(x, y2, label = "tangent line")
plt.plot(x, y, label = "line")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.show()
``````

# Python梯度圖

Python:

``````# 一次算x.y兩個偏微分(梯度)
h = 1e-4 # 0.0001
print(x.size)
for idx in range(x.size):

tmp_x = x[idx]

x[idx] = tmp_x + h
fxh1 = f(x)

x[idx] = tmp_x - h
fxh2 = f(x)

grad[idx] = (fxh1 - fxh2) / (2 * h)
x[idx] = tmp_x

# 計算所有x.y梯度

for idex, X in enumerate(x):

if __name__ == '__main__':
x = np.arange(-2.0, 2.0, 1.0)
y = np.arange(-2.0, 2.0, 3.0)

#將x和y擴展當相同大小.再轉一維
X, Y = np.meshgrid(x, y)
X = X.flatten()
Y = Y.flatten()

#求梯度

print(X)
print(Y)

#畫出quiver方向圖, xlim指定顯示大小, label標籤顯示, grid顯示網格
plt.xlim([-4, 4])
plt.ylim([-4, 4])
plt.xlabel('x0')
plt.ylabel('x1')
plt.grid()
plt.draw()
plt.show()
``````

x = [-2. -1. 0. 1. -2. -1. 0. 1.]
y = [-2. -2. -2. -2. 1. 1. 1. 1.]
x偏 = [-4. -2. 0. 2. -4. -2. 0. 2.]
y偏 = [-4. -4. -4. -4. 2. 2. 2. 2.]

# C#梯度圖

### 座標Class

``````class PointF2D
{
private float _x;
private float _y;

public PointF2D(float x, float y)
{
_x = x;
_y = y;
}

public PointF2D()
{
_x = 0;
_y = 0;
}

public float X
{
get
{
return _x;
}
set
{
_x = value;
}
}

public float Y
{
get
{
return _y;
}
set
{
_y = value;
}
}
}
``````

### 使用函數

``````interface Function
{
float Formula(float x, float y);
}

class Function1 : Function
{
public float Formula(float x, float y)
{
return x * x + y * y;
}
}
``````

### 微分函數

``````PointF2D Fun(Function fun, float x, float y)
{
float fun1 = 0.0f;
float fun2 = 0.0f;
float h = 1e-4f;

fun1 = fun.Formula((x + h), y);
fun2 = fun.Formula((x - h), y);
grad.X = (fun1 - fun2) / (h * 2);

fun1 = fun.Formula(x, (y + h));
fun2 = fun.Formula(x, (y - h));
grad.Y = (fun1 - fun2) / (h * 2);

}
``````

### 繪圖Class

1.因我們得知x向量y向量所以可以使用atan取得弧度(得到我們要畫的角度)。
2.之後我們可以先計算Cos和Sin用來取的邊長資訊。
3.取的x.y的方向(Direction)也就是要畫哪個象限。
4.maxLen計算要畫多長。
5.index是目前我們要畫的長度因此斜邊(bevel) = index / Cos。(Cos = 鄰/斜 推算)。
6.x乘上方向.y = Sin * 斜邊(Sin = 對/斜 推算)在乘上方向。
7.直到畫到我們目前最大長度。

``````class DrawF
{
PointF2D _center;
float _offset;

public DrawF(float offset)
{
_center = new PointF2D(5.0f * offset, 5.0f * offset);
_offset = offset;
}

public DrawF(float offset, float x, float y)
{
_center = new PointF2D(x * offset, y * offset);
_offset = offset;
}

public PointF2D getBlockPoint(float x, float y)
{
PointF2D point = new PointF2D();
point.X = (_center.X / _offset + x) * _offset;
point.Y = (_center.Y / _offset - y) * _offset;
return point;
}

public PointF2D Center
{
get
{
return _center;
}
}

public void drawBlock(Graphics graphics
, float x, float y)
{
Pen pen = new Pen(Color.FromArgb(255, 0, 0, 0), 1);
float xLen = x * _offset;
float yLen = y * _offset;

// 水平的分割出垂直數量
for (int row = 0; row <= x; row++)
{
graphics.DrawLine(pen, 0, row * _offset, yLen, row * _offset);
}

// 垂直的分割出水平數量
for (int col = 0; col <= y; col++)
{
graphics.DrawLine(pen, col * _offset, 0, col * _offset, xLen);
}
}

public void drawLine(Graphics graphics
, float xS, float yS
, float xE, float yE)
{
Pen pen = new Pen(Color.FromArgb(255, 178, 34, 34), 5);
graphics.DrawLine(pen, xS, yS, xE, yE);
}

public void drawPoint(Graphics graphics, Brush color
, float x, float y)
{
graphics.FillRectangle(color, x - 5.0f, y - 5.0f, 10.0f, 10.0f);
}

public void drawPointLine(Graphics graphics, Brush color
, float x, float y)
{
graphics.FillRectangle(color, x, y, 2.0f, 2.0f);
}

// point: 座標
// U: x向量
// V: y向量
// maxSum: 最大向量和
, PointF2D point
, float U, float V, float maxSum)
{
// 鄰邊 = U = x
// 對邊 = V = y
float absU = Math.Abs(U);
float absV = Math.Abs(V);

// atan取弧度(tan = 對邊 / 鄰邊
float radian = absU != 0 ? (float)Math.Atan(absV / absU) : (float)(Math.PI * 0.5);

// cos = 鄰邊 / 斜邊, sin = 對邊 / 斜邊
float xCos = U != 0.0f ? (float)Math.Cos(radian) : 1.0f;
float ySin = V != 0.0f ? (float)Math.Sin(radian) : 0.0f;

float xDirection = U < 0.0f ? -1.0f : U > 0 ? 1.0f : 0.0f;
float yDirection = V < 0.0f ? 1.0f : V > 0 ? -1.0f : 0.0f;

// 計算顯示長度比例
float max = absU > absV ? absU : absV;
float maxLen = (max / maxSum * _offset);

for (float index = 0; index < maxLen; index += 0.01f)
{
// 取得斜邊
// 取得x + 方向和y + 方向(對邊)
float bevel = index / xCos;
float x = index * xDirection;
float y = ySin * bevel * yDirection;

if (Math.Abs(y) > maxLen)
{
return;
}

drawPointLine(graphics, Brushes.Black, point.X + x, point.Y + y);
}
}

// point: 座標
// U: x向量
// V: y向量
, PointF2D point
, float U, float V)
{
float sum = Math.Abs(V) + Math.Abs(U);
}

// xyPoints: 座標陣列
// uvPoints: 向量陣列
, ArrayList xyPoints
, ArrayList uvPoints)
{
float maxUV = getMaxUV(uvPoints);

for (int index = 0; index < xyPoints.Count; index++)
{
PointF2D xyPoint = (PointF2D)xyPoints[index];
PointF2D uvPoint = (PointF2D)uvPoints[index];

}
}

// getMaxUV取得目前最大向量
// uvPoints: 向量陣列
private float getMaxUV(ArrayList uvPoints)
{
float maxUV = 0.0f;

for (int index = 0; index < uvPoints.Count; index++)
{
PointF2D uvPoint = (PointF2D)uvPoints[index];
float sum = Math.Abs(uvPoint.X) + Math.Abs(uvPoint.Y);
if (maxUV < sum)
{
maxUV = sum;
}
}

return maxUV;
}
}
``````

### Picture畫出

``````private void pictureBox1_Paint(object sender, PaintEventArgs e)
{
DrawF draw = new DrawF(50.0f, 5.0f, 5.0f);
draw.drawBlock(e.Graphics, 10.0f, 10.0f);
draw.drawPoint(e.Graphics, Brushes.Red, draw.Center);

ArrayList xyPoints = new ArrayList();
ArrayList uvPoints = new ArrayList();

float[] x = { -2, -1, 0, 1, -2, -1, 0, 1 };
float[] y = { -2, -2, -2, -2, 1, 1, 1, 1 };

Function1 fun = new Function1();

for (int index = 0; index < 8; index++)
{
PointF2D xyPoint = draw.getBlockPoint(x[index], y[index]);
PointF2D uvPoint = Fun(fun, x[index], y[index]);
Console.Write(uvPoint.X + " " + uvPoint.Y);
draw.drawPoint(e.Graphics, Brushes.Blue, xyPoint);
}