🏈指点迷津 | Brief🍇Python和MATLAB符号微分
计算导数的常用方法有 3 种:
数值微分依赖于导数的定义:
f′(x)=h→0limhf(x+h)−f(x) 其中,你把一个非常小的 h 和评估函数放在两个地方。这是最基本的公式,在实践中,人们使用其他公式,这些公式给出的估计误差更小。这种计算导数的方法主要适用于你不知道你的函数并且只能对其进行采样的情况。此外,对于高维函数,它需要大量的计算。
符号微分操纵数学表达式。如果你曾经使用过 matlab 或 mathematica,那么你会看到类似这样的内容
dxd(sin(x)x2cos(x−7))=x2sin(7−x)csc(x)+x2(−cos(7−x))cot(x)csc(x)+2xcos(7−x)csc(x) 在这里,他们知道每个数学表达式的导数,并使用各种规则(乘积规则、链式规则)来计算最终的导数。然后他们简化表达式以获得最终表达式。
自动微分可操纵计算机程序块。微分器具有对程序的每个元素求导的规则。它还使用链式法则将复杂表达式分解为更简单的表达式。
Python符号微分
导数规则要求一个项有 2 个子项,即二元运算(例如 *),以及一个一元运算(例如 sin)。
term.py:
from functools import reduce
class Term:
def __init__(self, label: str or int, children=None):
self.label = label
self.children = children if children else []
def __repr__(self):
if self.children:
return f"({self.label} " + reduce(lambda str1, str2: str1+" "+str2, (repr(child) for child in self.children)) + ")"
else:
return f"{self.label}"
@classmethod
def copy(cls, obj):
"""makes a copy of obj"""
return cls(obj.label, [cls.copy(child) for child in obj.children])
symbols_table:
MUL = '*'
ADD = '+'
SUB = '-'
DIV = '/'
POW = '^'
LOG = 'log'
SIN = 'sin'
COS = 'cos'
TAN = 'tan'
OPEN_BR = '('
CLOSE_BR = ')'
VAR = 'x'
symbols = {LOG, SIN, COS, TAN, MUL, ADD, SUB, DIV, POW, OPEN_BR, CLOSE_BR, VAR}
from term import Term
from symbols_table import *
import numbers
def calculate_derivative(term1: Term, rules: dict) -> Term:
if isinstance(term1.label, numbers.Number):
return Term(0)
elif term1.label == VAR:
return Term(1)
else:
return rules[term1.label](term1, rules)
def mul_derivative(term1: Term, rules) -> Term:
child1, child2 = term1.children
return Term(ADD, [Term(MUL, [calculate_derivative(child1, rules), child2]), Term(MUL, [child1, calculate_derivative(child2, rules)])])
def div_derivative(term1: Term, rules) -> Term:
child1, child2 = term1.children
return Term(DIV, [Term(SUB, [Term(MUL, [calculate_derivative(child1, rules), child2]), Term(MUL, [child1, calculate_derivative(child2, rules)])]),
Term(POW, [child2, Term(2)])])
def sum_derivative(term1: Term, rules) -> Term:
child1, child2 = term1.children
return Term(term1.label, [calculate_derivative(child1, rules), calculate_derivative(child2, rules)])
def power_derivative(term1: Term, rules) -> Term:
child1, child2 = term1.children
return Term(ADD,
[Term(MUL, [calculate_derivative(child2, rules), Term(LOG, [child1]), Term(POW, [child1, child2])]),
Term(MUL, [calculate_derivative(child1, rules), child2, Term(POW, [child1, Term(SUB, [child2, Term(1)])])])])
def sin_derivative(term1: Term, rules) -> Term:
return Term(MUL, [Term(COS, [term1.children[0]]), calculate_derivative(term1.children[0], rules)])
def cos_derivative(term1: Term, rules) -> Term:
return Term(MUL, [Term(-1), Term(SIN, [term1.children[0]]), calculate_derivative(term1.children[0], rules)])
def tan_derivative(term1: Term, rules) -> Term:
arg = Term.copy(term1.children[0])
return Term(ADD, [calculate_derivative(arg, rules), Term(MUL, [calculate_derivative(arg, rules), Term(POW, [term1, Term(2)])])])
def log_derivative(term1: Term, rules) -> Term:
return Term(DIV, [calculate_derivative(term1.children[0], rules), term1.children[0]])
rules_for_differentiation = {MUL: mul_derivative, ADD: sum_derivative, SUB: sum_derivative,
COS: cos_derivative, SIN: sin_derivative, TAN: tan_derivative,
POW: power_derivative, DIV: div_derivative, LOG: log_derivative}
def derivative(term1):
return calculate_derivative(term1, rules_for_differentiation)
MATLAB符号微分
求函数 y=f(x)^2 \frac{d f}{d x} 关于 f(x) 的导数
syms f(x) y
y = f(x)^2*diff(f(x),x);
Dy = diff(y,f(x))
\begin{aligned} & Dy = \\ & \qquad 2 f(x) \frac{\partial}{\partial x} f(x)\end{aligned}
找到描述质量弹簧系统运动的欧拉-拉格朗日方程。定义系统的动能和势能。
syms x(t) m k
T = m/2*diff(x(t),t)^2;
V = k/2*x(t)^2;
定义拉格朗日量。
L=2m(∂t∂x(t))2−2kx(t)2
欧拉-拉格朗日方程如下
0=dtd∂x˙∂L(t,x,x˙)−∂x∂L(t,x,x˙)
计算项 ∂L/∂x˙
D1 = diff(L,diff(x(t),t))
D1=m∂t∂x(t)
计算第二项∂L/∂x
D2(t)=−kx(t)
找到质量弹簧系统的欧拉-拉格朗日运动方程。
diff(D1,t) - D2 == 0
ans(t)=m∂t2∂2x(t)+kx(t)=0