第1题.同意楼上
第2题.修改后的代码,不知是不是你所要的结果:
clear
format long
syms k1 k01 VV % z=dφ/dt
V=[];
z=[];
for j=1:1:40
VVstep=0.0002; %VV代表电极电势φ
VVfinal=0.0160;
VV=0.008+VVstep*j
k1=150*exp(19.462*VV);
k01=150*exp(-19.462*VV);
equation_a=k01*200;
equation_b=-(5*k1+1000+k01*200);
equation_c=2*5*k1+1000+50*k01;
equation_d=-5*k1;
Coefficient=[equation_a equation_b equation_c equation_d]
temp=roots(Coefficient)
equation_y01=temp(1)
equation_x01=(200*equation_y01^2-200*equation_y01+50)/5
equation_z01=(0.059-VV)/(1.1*10^(-7))-96485*(k1*(1-equation_y01)^2-k01*equation_x01*equation_y01)/(10*10^(-10))
equation_y02=temp(2)
equation_x02=(200*equation_y02^2-200*equation_y02+50)/5
equation_z02=(0.059-VV)/(1.1*10^(-7))-96485*(k1*(1-equation_y02)^2-k01*equation_x02*equation_y02)/(10*10^(-10))
equation_y03=temp(3)
equation_x03=(200*equation_y03^2-200*equation_y03+50)/5
equation_z03=(0.059-VV)/(1.1*10^(-7))-96485*(k1*(1-equation_y03)^2-k01*equation_x03*equation_y03)/(10*10^(-10))
hold on
V=[V,VV];
z=[z;equation_z01,equation_z02,equation_z03];
end
plot(V,z)
legend('equation z01','equation z02','equation z03')
后注:
1.代码效率比较低,不容易阅读,变量名太长了点儿;
2.计算后有虚部产生,不知为何,能力有限。请检查算法是否无误
3.记得加分号以抑制中间结果的输出,如果想要输出中间结果最好用disp说明一下,以免阅读混乱;,