13518219792
Sympy是python中非常强大的符号运算库，可以以书写习惯表示数学表达式。下面介绍用Sympy求方程数值解的方法。
成都创新互联公司于2013年开始，是专业互联网技术服务公司，拥有项目成都做网站、网站建设网站策划，项目实施与项目整合能力。我们以让每一个梦想脱颖而出为使命，1280元永修做网站,已为上家服务,为永修各地企业和个人服务,联系电话:028-86922220
下面代码全部在
from sympy import *
init_printing(use_unicode=True) # 按书写习惯输出
下运行。
数学表达式的输入
首先声明符号：
x = symbols('x')
即计算机中的变量x代表数学表达式中的x。在后文输出中所有的x会显示为x。如果x=symbols('x0')，则输入的方程中所有x将在输出中以x0表示。
如果需要希腊字母
l, r = symbol('lambda rho')
l, r将分别以λ,ρ表示。可以在一个表达式中同时声明多个符号。
或者使用var()声明：
var('x')
与上面等效。
声明表达式：
f = (5/x)*(exp(x)-1)-exp(x)
此时若输出f可以看到书写习惯的表达式。由于表达式在markdown下显示不正常，在此不放置示例。注意f的类型是class 'sympy.core.add.Add'
求f(x)=0数值解
因为有的函数零点不止一个，因此在Sympy中解的输出为一个list。使用solve(表达式，自变量符号)可以解析地解方程：
s, = solve(f, x)
这里根据上面f的赋值，得到s为
LambertW(-5e**-5)+5
其中用了特殊函数表达。
我们需要求这个结果的数值近似，则输出
s.evalf()
得到输出
4.96511423174428
就是方程f(x)=0的数值解。
求给定自变量x值时函数f(x)的值 | 将表达式转化为函数
f.evalf(subs = {x:4.96})
得到f(4.96)的数值
0.141885450782171
如果需要以计算机函数的形式定义函数f(x)，则可以使用lambdify()进行转化：
f_func = lambdify(x, f)
之后可以调用
f_func(4.96)
输出
0.141885450782
利用这个方法可以测试方程的数值算法，如使用sympy接口写牛顿法等。
利用 numpy 很简单。可以利用pip安装
pip install numpy
然后（以你的方程为例），python 下
Python 2.7.10 (default, Oct 23 2015, 19:19:21)
[GCC 4.2.1 Compatible Apple LLVM 7.0.0 (clang-700.0.59.5)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
import numpy as np
a = np.array([[1,1],[1,-1]])
b = np.array([3,1])
print np.linalg.solve(a,b)
[ 2. 1.]
如果你学过 线性代数，那么这段代码很好理解。
不是很明确你需要做到什么程度，但基本可以通过以下两个手段得到：
手工解方程得到解析解，然后套入公式
使用一些工具包例如numpy可以自动求解
以下都给出例子
import numpy as np
import matplotlib.pyplot as plt
plt.axis("equal")
a = np.linspace(1,10,100) # a 的变化范围可以自己挑，前两个参数控制，
# 使用 numpy 自动求解
res = []
for x in a:
A = np.mat("1, 2; {}, -1".format(x))
b = np.mat("{}, 10".format(x)).T
res.append(np.linalg.solve(A, b))
# 计算完毕后取出每对x和y
x1 = [float(r[0]) for r in res]
y1 = [float(r[1]) for r in res]
plt.plot(x1, y1)
#####################################
# 手工计算过程很简单不放上来了，直接上结果
x2 = [(a1 + 20) / (2*a1 + 1) for a1 in a]
y2 = [(a1**2 - 10) / (2*a1 + 1) for a1 in a]
plt.plot(x2, y2)
本文归纳常见常微分方程的解析解解法以及基于Python的微分方程数值解算例实现。
考虑常微分方程的解析解法，我们一般可以将其归纳为如下几类：
这类微分方程可以变形成如下形式：
两边同时积分即可解出函数，难点主要在于不定积分，是最简单的微分方程。
某些方程看似不可分离变量，但是经过换元之后，其实还是可分离变量的，不要被这种方程迷惑。
形如
的方程叫做一阶线性微分方程，若 为0，则方程齐次，否则称为非齐次。
解法： (直接套公式)
伯努利方程
形如
的方程称为伯努利方程，这种方程可以通过以下步骤化为一阶线性微分方程：
令 , 方程两边同时乘以 ，得到
即 .
这就将伯努利方程归结为可以套公式的一阶线性微分方程。
形如
的方程称为二阶常系数微分方程，若 ，则方程称为齐次的，反之称为非齐次的。以下默认方程是非齐次的。
求解此类方程分两步：
原方程的解 = 齐次通解 + 非齐次特解
首先假设 .用特征方程法，写出对应的特征方程并且求解：
解的情况分为以下三种：
情况一：方程有两个不同的实数解
假设两个实数解分别是 , 此时方程的通解是
情况二：方程有一个二重解
假设该解等于 ，此时方程的通解是
情况三：方程有一对共轭复解
假设这对解是 , 此时方程的通解是
对于 和特征根的情况，对特解的情况做如下归纳：
形如
的方程叫做高阶常系数微分方程，若 ，则方程是齐次的，否则是非齐次的。下面默认方程是非齐次的。
求解此类方程分两步：
原方程的解 = 齐次通解 + 非齐次特解
考虑带有第三类边界条件的二阶常系数微分方程边值问题
问题一:两点边值问题的解析解
由于此方程是非齐次的,故 求解此类方程分两步：
原方程的解 = 齐次通解 + 非齐次特解
首先假设 . 用特征方程法，写出对应的特征方程
求解得到两个不同的实数特征根: .
此时方程的齐次通解是
由于 . 所以非齐次特解形式为
将上式代入控制方程有
求解得: , 即非齐次特解为 .
原方程的解 = 齐次通解 + 非齐次特解
于是,原方程的全解为
因为该问题给出的是第三类边界条件,故需要求解的导函数
且有
将以上各式代入边界条件
解此方程组可得: .
综上所述,原两点边值问题的解为
对一般的二阶微分方程边值问题
假定其解存在唯一,
为求解的近似值, 类似于前面的做法,
考虑带有第三类边界条件的二阶常系数微分方程边值问题
问题二:有限差分方法算出其数值解及误差
对于 原问题 ,取步长 h=0.2 ,用 有限差分 求其 近似解 ,并将结果与 精确解y(x)=-x-1 进行比较.
因为
先以将区间划分为5份为例，求出数值解
结果:
是不是解出数值解就完事了呢？当然不是。我们可以将问题的差分格式解与问题的真解进行比较，以得到解的可靠性。通过数学计算我们得到问题的真解为 ，现用范数来衡量误差的大小：
结果:
接下来绘图比较 时数值解与真解的差距：
结果:
将区间划分为 份, 即 时.
结果:
绘图比较 时数值解与真解的差距：
最后，我们还可以从数学的角度分析所采用的差分格式的一些性质。因为差分格式的误差为 ， 所以理论上来说网格每加密一倍，与真解的误差大致会缩小到原来的 . 下讨论网格加密时的变化：
结果：
方法/步骤
用Python解数学方程，需要用到Python的一个库——SymPy库。
SymPy是符号数学的Python库，它的目标是成为一个全功能的计算机代数系统，同时保持代码简洁、易于理解和扩展。
如果你的电脑上还没有安装sympy库，那就赶紧安装吧，安装命令：
pip3 install sympy
请点击输入图片描述
先来解一个简单点的方程吧。
题目： 5x + 20 = 100
先直接上代码：
from sympy import *
x = Symbol('x')
print(solve([5*x + 20 - 100], [x]))
请点击输入图片描述
再来一个复杂点的二元一次方程吧。
题目：3x + 4y =49, 8x- y = 14
代码如下：
from sympy import *
x = Symbol('x')
y = Symbol('y')
print(solve([3*x + 4*y - 49, 8*x - y - 14], [x, y]))
请点击输入图片描述
有没有发现规律呢，简单总结一下：
1）变量赋值，使用symbol函数转换；
2）将方程式移到方程的左边，使右边等于0；
3）使用solve函数解方程。
当然了，python的基础语法必须掌握，至少需要掌握python最基础的算数运算符。
+ 加 ---- 两个对象相加
- 减 ----- 得到负数或是一个数减去另一个数
* 乘 ----- 两个数相乘或是返回一个被重复若干次的字符串
/ 除 ----- x 除以 y
% 取模 ----- 返回除法的余数
** 幂 ----- 返回x的y次幂
log() 对数-----对数 log()
下面来个难度大点的方程。
请点击输入图片描述
代码如下：
from sympy import *
t = Symbol('t')
x = Symbol('x')
m = integrate(sin(t)/(pi-t), (t, 0, x))
print(integrate(m, (x, 0, pi)))
请点击输入图片描述
