Matlab验证：谭编首次提出的曲线微分求导方法科学吗？
一劳永逸地搞定微分电容
Crusher
武汉大学博士生
能源材料是一个很热的话题，而对于奋战在一线的研究者们来说，循环伏安（CV）和线性伏安（LSV）测试几乎必不可少。
然而，在大量文献中，CV曲线似乎并非必需的数据曲线，取而代之的是微分电容曲线。
不少同学参考这些文献方法，试图得到类似的微分容量曲线，但是，实际上会遇到很多困难。
要么不知道怎么得到（当然某些国产仪器可以从充放电曲线中点击查看微分容量曲线，某些进口电化学工作站未必会有这个功能），要么在利用Origin绘图软件进行微分处理时，得到的一定是类似下图的奇怪结果。
图1 采用Origin微分得到的微分容量曲线
满屏的竖线，第一感觉往往是：头大！
至于为什么会这样，主要是由于采点的原因，或者仪器精密度的限制，导致直接作微商的过程出现了很多0/0的不定型，而0/0的未定型在数学上即可以是0，也可以是无穷大或无穷小。
今天谭编联合武汉大学的某博士科研狗，给大家带来了福音：一劳永逸地搞定微分电容。
在今天发布的【微分求导】专题系列文章中，小编由主到次、由浅入深解决以下问题:
问题1：怎样求微分电容电容？（本文）
问题2：微分电容到底是什么？（下一篇）
问题3：从数学上进一步深层探讨求解微分电容方法的可信性和局限性。（下一篇）
问题1，怎样求微分电容
本文的具体验证原理基于【编辑之谭】前期发布的《Origin+Excel神操作：从充放电曲线绘制微分电容曲线》，请点击以下链接查阅。
Origin+...神操作：从充放电曲线绘制微分电容曲线
在上面这篇文章中，谭编首次提出了微分电容的平均求导法（方法1）、交叉求导法（方法2）。
本文对这两种方法的科学性进行验证，同时，综合考虑方法1和方法2有可能导致正负偏差或计算结果有可能引入的曲线极化现象，本文提出均值法（方法3），即对两种方法的计算结果再次求均值的方法，对比研究并验证这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))
10holdon
11end
12forj=1:n/4
13name(2*j-1,:)=['第'num2str(j),'次充电'];
14name(2*j,:)=['第'num2str(j),'次放电'];
15legend(name)
16title('充放电曲线')
17end
18d=150;
19figure(2)
20fori=1:length(A(:,1))-d
21forj=1:n/2
22C(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)));
24end
25end
26forj=1:n/2
27plot(C(:,2*j-1),C(:,2*j))
28holdon
29end
30forj=1:n/4
31name1(2*j-1,:)=['第'num2str(j),'次充电微分电容'];
32name1(2*j,:)=['第'num2str(j),'次放电微分电容'];
33legend(name1)
34title('方法1：平均求导法')
35end
36d1=150;
37figure(3)
38fori=1:length(A(:,1))-d1
39forj=1:n/2
40D(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)));
42end
43end
44forj=1:n/2
45plot(D(:,2*j-1),D(:,2*j))
46holdon
47end
48forj=1:n/4
49name2(2*j-1,:)=['第'num2str(j),'次充电微分电容'];
50name2(2*j,:)=['第'num2str(j),'次放电微分电容'];
51legend(name2)
52title('方法2：交叉求导法')
53end
54%方法3将两者相加
55figure(4)
56fori=1:min(length(A(:,1))-d,length(A(:,1))-d1)
57forj=1:n/2
58E(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));
60end
61end
62forj=1:n/2
63plot(E(:,2*j-1),E(:,2*j))
64holdon
65end
66forj=1:n/4
67name3(2*j-1,:)=['第'num2str(j),'次放电微分电容'];
68name3(2*j,:)=['第'num2str(j),'次充电微分电容'];
69legend(name3)
70title('方法3：两者方法取均值')
71end
72figure(5)
73forj=1:n/2%同时比较三种方法求出的第一周充放电曲线的微分电容
74plot(C(:,2*j-1),C(:,2*j),'g')
75holdon
76plot(D(:,2*j-1),D(:,2*j),'b')
77holdon
78plot(E(:,2*j-1),E(:,2*j),'m')
79legend('方法1','方法2','方法3')
80title('三种方法比较')
81end
不需要懂电化学知识和数学原理，尤其矩阵运算方面的知识，只需要将你的数据拷入到excel表格中，按电压-容量（或比容量，都可以二者无所谓），有一条曲线就粘贴一次“电压-容量”，有多条曲线就粘贴多次。多说无益，我们直接上图吧。
图2 将充放电数据拷入Excel文件
为了验证脚本的可信性，我们分峰形比较简单和多对峰分别讨论。
1.简单峰（LiFePO4）
LiFePO4是一个成熟的电极材料，目前被比亚迪，宁德时代等各大电池厂商大量使用。磷酸铁锂有一个平坦的工作平台，但可逆容量不甚高，下图是我们曾经做的一次掺杂，可见充放电平台基本上保持不变。现在利用我们编写的Matlab脚本程序，自动读取该曲线的Excel数据文件，计算它的微分电容曲线。
图3 改性的磷酸铁锂的充放电曲线(研究对象)
利用我们编写的Matlab脚本，对方法1、方法2、方法3进行了验证，其结果如图4所示。
图4 三种方法的微分容量曲线
（点击图片，可放大查阅）
可见，三种方法几乎没有明显差别（放电的小峰正是由于掺杂引起的），这主要是步长d，d1都取的较小。如果得不到较好的结果，请直接尝试修改不同的步长值。另外，如果步长较大时，推荐使用方法3，因为方法1和2会分别产生一个负的和正的误差，两者叠加能在一定程度上消除数据处理带来的误差。
上述的验证结果表明：我们编写的Matlab脚本对简单峰的识别是非常成功的，那么如果是多条曲线，而且是多个充放电平台（多峰）曲线的复杂情况呢？
2. 多峰多曲线
图5是一种具有多个充放电平台的充放电曲线，不用改脚本，直接将充放电的数据拷入到excel表格中，用matlab读入后，运行即可。
运行结果：
图5 具有多平台的充放电曲线（研究对象）
图6 多平台充放电曲线的微分求导结果
（点击图片，可放大查看）
对比可见，除了主要的峰形，3.2V等处的小峰也清晰可见，这有助于对充放电机理深刻理解。所以该脚本对于多对峰，多条曲线的处理能力也是很不错的。唯一要注意的是，如果d取的较大，可能会造成曲线的左右移动，人为地制造出极化，为此，可将第39行（即第5行，下面只是节选代码）的d1除以2。
1d1=150;
2figure(3)
3fori=1:length(A(:,1))-d1
4forj=1:n/2
5D(i,2*j-1)=A(i+d1/2,2*j-1);
6D(i,2*j)=((A(i+d1,2*j))-(A(i,2*j)))/((A(i+d1,2*j-1))-(A(i,2*j-1)));
7end
8end
说了这么多，如果某些同学习惯在origin中作图，怎样导出微分后的数据呢，只需加一条输出语句‘xlswrite’即可，具体用法请找“度娘”（百度）。
如果只是在文章中单纯地用一下，上面的内容就够了。但为了深刻理解，谭编推荐阅读本期专题的其他文章。