matlab 1到无穷_Matlab：一劳永逸搞定微分电容

4月25日，『编辑之谭』分享了两种微分求导方法(平均求导方法、交叉求导方法)，一出就引起了广泛关注。

Origin ...神操作：从充放电曲线绘制微分电容曲线

(感谢Crusher指出笔误：斜率应为tan())

谭编首次提出的两种微分求导方法可以处理电池的充放电曲线数据，计算微分电容曲线数据，完美绘制电池的微分容量曲线，弥补Origin绘图软件微分不成功的缺陷。

那么这两种微分方法到底呢？

科不科学？严不严谨？可不可行？

有没有修正错误的方法？

看看武汉大学Crusher博士怎么用Matlab编程来验证这些问题！

(满足各大需求学术公众号需要转载，将上一期专题中的系列文章汇编成本文，欢迎联系谭编转载。)

一劳永逸地完成微分电容

Crusher

武汉大学博士生

能源材料是一个热门话题，对于奋战在一线的研究人员来说，伏安的循环(CV)和线性伏安(LSV)测试几乎是必不可少的。然而，在大量的文献中，CV曲线似乎不是必要的数据曲线，而是微分电容曲线。

许多学生参考这些文献方法，试图获得类似的微容量曲线，但实际上会遇到很多困难。要么不知道如何得到它，要么使用它Origin微分处理绘图软件时，必须得到类似下图的奇怪结果。

图1 采用Origin微分得到的微分容量曲线

满屏竖线，第一感觉往往是:头大！

至于为什么会这样，主要是因为仪器采集点或者仪器精度的限制，导致直接作为微信业务的过程中出现了很多0/0的不定型，而0/0的不定型可以是数学上的0，也可以是无限小。

虽然谭编提出的两种方法完美地获得了微分电容曲线，但这两种方法的科学性质引起了我们的质疑。作者采用Matlab编程，从Excel在数据文件中读取充放电数据，一键处理并绘制微分电容曲线。

谭编提出的平均导数（方法1）和交叉导数（方法2），以及作者综合考虑这两种方法的误差提出的平均值法（方法3），即对两种方法的计算结果进行比较研究，验证电池充放电曲线的科学性、准确性和可行性。同时，一键处理磷酸铁锂单平台、多平台特征的充放电数据和多个循环数据，验证其可行性。

本文主要讨论以下3个问题：

问题1:如何求微分电容？

问题2:什么是微分电容？

问题3：进一步探讨微分电容解决方案的可信度和局限性。

如何求微分电容？

我们编写的Matlab脚本如下(左右滑动，可阅读完整代码)，可复制Matlab源代码，也可留言或联系谭编微信免费获取。


1%Matlab源代码 2clc 3clearall 4A=xlsread('testnfpcv.xlsx'); 5n=size(A,2) 6A; 7figure(1) 8fori=1:n/2 9plot(A(:,2*i),A(:,2*i-1))10holdon11end12forj=1:n/413name(2*j-1,:)=['第'num2str(j),14.name(2*j,:)=['第'num2str(j),二次放电；15.legend(name)16title17.充放电曲线end18d=150;19figure(2)20fori=1:length(A(:,1))-d21forj=1:n/222C(i,2*j-1)=(A(i,2*j-1) A(i d,2*j-1))/2;23C(i,2*j)=((A(i d,2*j))-(A(i,2*j)))/((A(i d,2*j-1))-(A(i,2*j-1)));24end25end26forj=1:n/227plot(C(:,2*j-1),C(:,2*j))28holdon29end30forj=1:n/431name1(2*j-1,:)=['第'num2str(j),三二次充电微分电容器name1(2*j,:)=['第'num2str(j),微分电容电微分电容legend(name1)34title35end36d1=150;37figure(3)38fori=1:length(A(:,1))-d139forj=1:n/240D(i,2*j-1)=A(i d1/2,2*j-1);41D(i,2*j)=((A(i d1,2*j))-(A(i,2*j)))/((A(i d1,2*j-1))-(A(i,2*j-1)));42end43end44forj=1:n/245plot(D(:,2*j-1),D(:,2*j))46holdon47end48forj=1:n/449name2(2*j-1,:)=['第'num2str(j),500次充电微分电容器name2(2*j,:)=['第'num2str(j),51legend(name2)52title53end54%方法3将两者相加55figure(4)56fori=1:min(length(A(:,1))-d,length(A(:,1))-d1)57forj=1:n/258E(i,2*j-1)=0.5*(C(i,2*j-1) D(i,2*j-1));59E(i,2*j)=0.5*(C(i,2*j) D(i,2*j));60end61end62forj=1:n/263plot(E(:,2*j-1),E(:,2*j))64holdon65end66forj=1:n/467name3(2*j-1,:)=['第'num2str(j),68name3(2*j,:)=['第'num2str(j),69legend(name3)70title711end72figure(5)73forj=1:n/2%同时，比较第一周充放电曲线的微分电容器74plot(C(:,2*j-1),C(:,2*j),'g')75holdon76plot(D(:,2*j-1),D(:,2*j),'b')77holdon78plot(E(:,2*j-1),E(:,2*j),'m')79legend方法1'，方法2'，'方法3')80title811end

不需要了解电化学知识和数学原理，特别是矩阵操作知识，只需复制您的数据excel在表中，根据电压-容量（或比容量，两者都不重要），有一条曲线粘贴一次电压-容量，有多条曲线粘贴多次。多说是没有用的，让我们直接在上面的图片中。

图2  将充放电数据拷入Excel文件

为了验证脚本的可信性，我们分峰形比较简单和多对峰分别讨论。

1.1 简单峰(LiFePO4) 

LiFePO4 是一种成熟的电极材料，目前被比亚迪，宁德时代等各大电池厂商大量使用。磷酸铁锂有一个平坦的工作电压平台，但可逆容量不甚高，下图是我们曾经做的一次掺杂，可见充放电平台基本上保持不变。现在利用我们编写的Matlab脚本程序，自动读取该曲线的Excel数据文件，计算它的微分电容曲线。

图3  改性的磷酸铁锂的充放电曲线(研究对象)

利用我们编写的Matlab脚本，对方法1、方法2、方法3进行了验证，其结果如图4所示。

图4  三种方法的微分容量曲线

(点击图片，可放大查阅)

可见，三种方法几乎没有明显差别(放电的小峰正是由于掺杂引起的)，这主要是步长d，d1都取得较小。如果得不到较好的结果，请直接尝试修改不同的步长值。另外，如果步长较大时，推荐使用方法3，因为方法1和2会分别产生一个负的或正的误差，方法3采用两者均值，能在一定程度上可以消除数据处理带来的误差。

(这是一个广告)

上述的验证结果表明：我们编写的Matlab脚本对简单峰的识别是非常成功的，那么如果是多条曲线，而且是多个充放电平台(多峰)曲线的复杂情况呢？

1.2 多峰多曲线

图5是一种具有多个充放电平台的充放电曲线，不用改脚本，直接将充放电的数据拷入到excel表格中，用Matlab读入后，运行即可。

运行结果：

图5  具有多平台的充放电曲线(研究对象)

图6  多平台充放电曲线的微分求导结果

(点击图片，可放大查看)

对比可见，除了主要的峰形，3.2V等处的小峰也清晰可见，这有助于对充放电机理的深刻理解。所以该脚本对于多对峰、多条曲线的处理能力也是很不错的。唯一要注意的是，如果d取得较大，可能会造成曲线中峰的起止位置发生左右漂移，人为地制造出极化，为此，可将第39行(即第5行，下面只是节选代码)的d1除以2。

1d1=150;2figure(3)3for i=1:length(A(:,1))-d14    for j=1:n/25    D(i,2*j-1)=A(i+d1/2,2*j-1);6    D(i,2*j)=((A(i+d1,2*j))-(A(i,2*j)))/((A(i+d1,2*j-1))-(A(i,2*j-1)));  7    end8end

说了这么多，如果某些同学习惯在origin中作图，怎样导出微分后的数据呢，只需加一条输出语句‘xlswrite’即可，具体用法请找“度娘”(百度)。

我们提供了这个Matlab脚本，欢迎在本文下方留言，看大家的反响，如果有必要的话，后续谭编再出一个Python脚本。

如果你是个实用主义青年，只是在文章中单纯地用一下，上面的内容就够了，学到这一步就可以了，这个Matlab脚本足以应付此类曲线的数据处理。

但勤奋好学的你肯定不满足于此，肯定想要理解微分电容底层的电化学知识，那么请往下看。

2 微分电容到底是什么？

细心地人可能早就发现，dQ/dV~V的峰形跟cv的是一样的，但具体为什么可能有些人并不理解，讲得稍微深刻一点的地方，一般是这么介绍的：

(公式采用【编辑之谈】公众号的【upub编辑器】)

后面呢？没了。那为什么微分电容和CV的结果二者一致呢？

事情得从一个原始的谁都知道，但绝大多数没怎么理解的东西说起，这个东西叫：平衡电势。

什么是平衡电势，百度上的定义是这样的：

这段话大家都觉得很有道理，也都读得很明白。那还是这个问题，(用人话讲)平衡电势到底是什么？

平衡电池指的是在某种状态下，一个物质(设为M)不从外回路得到/失去电子，变成自己的还原/氧化态(也就是M坚持做自己)的情况下，材料自由能的外在体现。

如果外部条件变化，这种平衡就会被打破，这正好就是贯通热力学和动力学的桥梁：

(公式采用【编辑之谈】公众号的【upub编辑器】)

为什么这么说呢？如果有外力，比如出现另外一个比它的氧化态还强的氧化剂，或比还原态还原性还强的还原剂，或者直接充放电(本质还是强行从材料夺电子或给电子)。

更形象点说，当外部氧化剂(设为O)的电位比M的高时，M就会被O夺走电子(热力学上说是这样，但具体能不能实现还要看这个过程有没有其他阻力，所以就有了所谓动力学的问题)，这时候是不做功的.

(这还是一个广告)

但如果将二者隔开(设法隔开)，并在中间架一根导线，就能让电子从外部流动，这就是成了电池。 

但如果是人为地缓慢的提高电势，而不是直接给一个强的氧化剂去夺取电子，会怎么样呢？

开始时，电位降低，电子夺取不下来，外电路电流很小(从而电量，即容量也很小，q=it)。

当电位提高到到平衡电势时并稍加一个过电位克服极化(这可能又是一个谁都知道但大多数人理解得一般的概念，以后有机会开专题讲)的阻力，这时候M就会快速地变成它的氧化态，剧烈地向外给出电子，放电过程则正好相反。

如果这个电子出得去、回得来，那就可以做成二次电池反复利用；如果出得去回不来，那就只能做成一次干电池。

有了上面的概念，我们再来理解为什么微分电容和CV是一回事呢？

上面提到，人为地提高电势到平衡电位以上，会剧烈地夺取M上的电子，对外显示的就是电流/容量。而认为提高电势，即

(1)以一个横定的速度直接提高电势(CV)；

(2)以一个横定的速度灌入电流(电势不是匀速上升)来实现，这就是恒电流充放电曲线。当然也可以既非恒电流也非恒电势，但不便于研究问题，用的不多。

如果是(1)，当电势低于平衡电势时，M不会受影响，不会失去自身的电子。当电势到达平衡电势时，这是由于外力的存在，M体系内的能量发生变化，电子被外回路夺取，对外显示，就是出现了峰。如果走1圈，那就是循环伏安CV，只走半圈，那就是线性伏安LSV。所以本质上说，CV中出现峰值，是说材料在这个电位处被打破了平衡。至于出现多个峰，即可能由于不同电对中心(如Fe和Mn)的平衡电势本就不一样，也可能由于同一个电对(比如Fe)在材料中处的晶体位置不一样，这个位置影响了它的平衡电势。[1]

如果是(2)，向M不断灌入电流时，电位也在发生变化，如果低于平衡电势，灌入的电流相当于直接给了一个导体正电荷，但它无法吸纳，就会使得的电位迅速变化，这就是充放电曲线中的类似于直线的部分(如图7中被蓝色框内的部分)。但当电势继续升高，接近第一个平衡电势，就会出现正电荷被吸收，导致电压升高得较慢，如第一个圈。如果第二个平衡电势处的M特别多，就会出现给多少正电荷都能吸收(M处于吃正电荷吃不饱的状态)，就会使得电势基本维持了一个平台，而这时候，还是在向外回路剧烈给电子。

图 7  充放电曲线的特征分析

到这里就可以看出：不管是CV还是重放电曲线，充电到达平衡电势并超过它时，都会向外输出大量电子，电子最直观的体现就是电流增加，所以在CV曲线中会出现一个峰值，而在充放电曲线中会出现一个平台，二者的本质是完全一致的。做一个简单的微分变化：

(公式采用【编辑之谈】公众号的【upub编辑器】)

平衡电势负极处，会伴随着一个很大的电流I。

还有一种理解方式是：将上面的图形逆时针旋转90°(如下)，即做成容量-电压曲线，你们发现了什么？ 正是Q-V 曲线在每个V 处的斜率，而平台处正是曲线斜率最大的地方！

图 8  充放电曲线上斜率的物理意义

到这里，小编将手把手地帮你彻底了解了微分电容是怎么一回事，从而做到库里有粮，心里不慌，以不变应万变。如果还有不理解的地方欢迎在本文末尾留言提问。

如果你是个标准的化学学霸而对数学不感兴趣，那么到这里也就可以了，但是作为祖国新时代的大学生，我想你会忍不住看第3个问题。

3 我们提出的曲线微分方法数学上严谨吗？

答案是：不严谨，但是……。

导数和微分直接的关系是什么？变化量和微分之间的关系又是什么？

我们先看导数定义：

其中要求Δx趋于0，而这在我们采的数据点中是不成立的。导数也可以这样定义：

图9  微分的定义及其几何意义

可见，导数是微分的比值。而微分与变化量之间的关系为：

图10  微分与变化量之间的关系

变化量Δy等于微分dy以及和一个比dy更高阶的无穷小的和，所以dy只是Δy的主要部分。而在我们的dQ/dV 曲线的求法中，我们直接以ΔQ代替了dQ，这是不严谨的。

同时还有另一个地方也有问题：因为避免出现竖线(即出现很多0/0的不定型)，我们对曲线上的点都做了一定的处理，比如跨步长取点。

这样实际上求出来的是两个点之间割线的斜率，而不是某一个点处的斜率。

那怎样完美的求微分电容曲线呢？理想的答案是有一台非常非常精密的充放电仪(目前一般实验室用的都在微安级别)，做出十分光滑的充放电曲线，直接对该曲线严格按数学方法求导，便是真正的微分电容曲线。

这可能是“你想”而不是“理想”的状态吧？！

尽管如此，该文还是为计算微分电容曲线提供了可行的方法，在一般的数据处理中，已经足够了。

【参考文献】

[1] X. Pu, H. Wang,T. Yuan, S. Cao, S. Liu, L. Xu, H. Yang, X. Ai, Z. Chen, Y. Cao, Energy Stor.Mater. (2019) DOI 10.1016/j.ensm.2019.1002.1017.

4 谭编怎么说？

本文提出的这3种微分求导方法，可以针对任何曲线，计算其微分曲线，可以从微分曲线中找到非常细微的变化。例如，从热重(TG)曲线数据，通过计算微分曲线得到相变温度、结晶水失水反应温度、脱水反应温度、分解反应温度、恒重温度等等。大家可以尝试验证一下，从TG曲线数据计算其微分曲线，并与实验测得的一次微分DTG曲线对比，研究是否吻合？如果是，那我们可以不用做DTG曲线了！(哈哈，想得美！)

有兴趣的网友，欢迎利用这两种方法，对您所遇到的任何实验数据曲线进行微分处理，验证一下谭编提出的这两种微分求导方法的可行性、普适性。

欢迎广大网友将你的验证结果投稿到『编辑之谭』，谭编自费为您发放百十来块的稿酬哦！

欢迎在本文末尾留言，若您留言提出的绘图问题收到的点赞最多，谭编将第一时间将该问题的解决教程发布出来。