以下是程序:
clear;
clc;
E=20
L=3.3e-3
C=1e-3
R=19;
T=4e-4;
a = 1/(2*R*C);
w = sqrt(1/L/C-1/4/R^2/C^2);
A2 = [0,-1/L;1/C,-1/R*C];
A1 = A2;
A3=[0,0;0,-1/C*R];
B1=[1/L;0];
B2=[0;0];
B3=B2;
I=[1,0;0,1];
for m=1:5001
Iref=0.1+0.01*(m-1);
i(1)=0;
u(1)=0;
for n=1:5000
X = [ i( n ); u( n ) ];
p1=(exp(-a*T)/w*((a*sin(w*T))+(w*cos(w*T))));
q1=exp(-a*T)/w*sin(w*T);
Ib1=(Iref/p1)+(q1*u(n)/L/p1)-(((1-p1)/R)+(q1/L))*E/q1;
if i(n)<=Ib1
X(:,n+1)=exp(A1*T)*X+(exp(A1*T)-I)*A1^-1*B1*E;
else
syms t1;
p3=exp(A1*t1)*X+(exp(A1*t1)-I)*A1^-1*B1*E;
q3=p3(1,:);