matlab length_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 3clear all 4A=xlsread('test nfpcv.xlsx'); 5n=size(A,2)  6A; 7figure(1) 8for i=1:n/2 9    plot(A(:,2*i),A(:,2*i-1))10    hold on11end12for j=1:n/4 13    name(2*j-1,:)=['第' num2str(j),'次充电'];14    name(2*j,:)=['第' num2str(j),'次放电'];15    legend(name) 16    title('充放电曲线')17end18d=150;19figure(2)20for i=1:length(A(:,1))-d21    for j=1:n/222    C(i,2*j-1)=(A(i,2*j-1)+A(i+d,2*j-1))/2;23    C(i,2*j)=((A(i+d,2*j))-(A(i,2*j)))/((A(i+d,2*j-1))-(A(i,2*j-1)));      24    end25end26for j=1:n/227    plot(C(:,2*j-1),C(:,2*j))28    hold on29end30for j=1:n/4 31    name1(2*j-1,:)=['第' num2str(j),'次充电微分电容'];32    name1(2*j,:)=['第' num2str(j),'次放电微分电容'];33    legend(name1) 34    title('方法1：平均求导法')35end36d1=150;37figure(3)38for i=1:length(A(:,1))-d139    for j=1:n/240    D(i,2*j-1)=A(i+d1/2,2*j-1);41    D(i,2*j)=((A(i+d1,2*j))-(A(i,2*j)))/((A(i+d1,2*j-1))-(A(i,2*j-1)));  42    end43end44for j=1:n/245    plot(D(:,2*j-1),D(:,2*j))46    hold on47end48for j=1:n/4  49    name2(2*j-1,:)=['第' num2str(j),'次充电微分电容'];50    name2(2*j,:)=['第' num2str(j),'次放电微分电容'];51    legend(name2) 52    title('方法2：交叉求导法')53end54%方法3 将两者相加 55 figure(4)56for i=1:min(length(A(:,1))-d,length(A(:,1))-d1)57    for j=1:n/258    E(i,2*j-1)=0.5*(C(i,2*j-1)+D(i,2*j-1));59    E(i,2*j)=0.5*(C(i,2*j)+D(i,2*j));   60    end61end62for j=1:n/263    plot(E(:,2*j-1),E(:,2*j))64    hold on65end66for j=1:n/4  67    name3(2*j-1,:)=['第' num2str(j),'次放电微分电容'];68    name3(2*j,:)=['第' num2str(j),'次充电微分电容'];69    legend(name3) 70    title('方法3：两者方法取均值')71end72figure(5)73for j=1:n/2 % 同时比较三种方法求出的第一周充放电曲线的微分电容74    plot(C(:,2*j-1),C(:,2*j),'g')75    hold on76    plot(D(:,2*j-1),D(:,2*j),'b')77    hold on78    plot(E(:,2*j-1),E(:,2*j),'m')79    legend('方法1','方法2','方法3')80    title('三种方法比较')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)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’即可，具体用法请找“度娘”(百度)。

如果只是在文章中单纯地用一下，上面的内容就够了。但为了深刻理解，谭编推荐阅读本期专题的其他文章。

谭编给大家推荐几本参考书，可以点击后去当当、京东等平台自行购买。

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

科学绘图，图绘科学！

如果本文对您的科学绘图有用，请记得转发朋友圈和微信群！

编辑之谭