laplace 偏微分方程 MATLAB help

clc,clear all
syms L C I1(t) Q(t) s
R = sym('R%d',[1 3])    %矩阵变量
assume([t L C R] > 0)  %positive
E(t) = 1*sin(t)        % Voltage = 1 V
微分方程
dI1=diff(I1,t)
dQ = diff(Q,t)
eqn1 = dI1+(R(2)/L)*dQ == (R(2)-R(1))/L*I1
eqn2 = dQ == (1/(R(2)+R(3))*(E-Q/C)) + R(2)/(R(2)+R(3))*I1

假设电流和电荷都是0
cond1 =I1(0)==0
cond2 = Q(0)==0

对方程进行laplace变换
eqn1LT =laplace(eqn1,t,s)  %t是需要进行变换的变量，s是默认的符号，可以改
eqn2LT =laplace(eqn2,t,s)  %t是需要进行变换的变量，s是默认的符号，可以改

用函数solve来求解符号变量
函数solve仅求解符号变量。因此, 要使用solve, 首先 substitute laplace(I1(t),t,s) and laplace(Q(t),t,s) with the variables I1_LT and Q_LT
syms I1_LT Q_LT
eqn1LT =subs(eqn1LT,[laplace(I1,t,s) laplace(Q,t,s)],[I1_LT Q_LT])
eqn2LT = subs(eqn2LT,[laplace(I1,t,s) laplace(Q,t,s)],[I1_LT Q_LT])
求解I1_LT和Q_LT的方程.
eqns = [eqn1LT eqn2LT];
vars = [I1_LT Q_LT];
[I1_LT, Q_LT] = solve(eqns,vars)
计算laplace的逆变换得到Q，并且化简结果
I1sol = ilaplace(I1_LT,s,t);
Qsol = ilaplace(Q_LT,s,t);
I1sol = simplify(I1sol);
Qsol = simplify(Qsol);

在绘制结果之前, 请用电路元素的数值替换符号变量
R1 = 4 Ω , R2 = 2 Ω, R3 = 3 Ω, C = 1/4 F, L = 1.6 H , I1 (0)= 15 A , and Q(0) = 2 C.
vars = [R L C I1(0) Q(0)];
values = [4 2 3 1.6 1/4 15 2];
I1sol = subs(I1sol,vars,values)
Qsol = subs(Qsol,vars,values)
plot results
Show both the transient and steady state behavior by using two different time intervals: 0 ≤ t ≤ 10 and 5 ≤ t ≤ 25
subplot(2,2,1)
fplot(I1sol,[0 10])                      
title('Current')
ylabel('I1(t)')
xlabel('t')

subplot(2,2,2)
fplot(Qsol,[0 10])                        
title('Charge')
ylabel('Q(t)')
xlabel('t')

subplot(2,2,3)
fplot(I1sol,[5 25])                  
title('Current')
ylabel('I1(t)')
xlabel('t')
text(7,0.25,'Transient')
text(16,0.125,'Steady State')

subplot(2,2,4)
fplot(Qsol,[5 25])                        
title('Charge')
ylabel('Q(t)')
xlabel('t')
text(7,0.25,'Transient')
text(15,0.16,'Steady State')


Initially, the current and charge decrease exponentially. However, over the long term, they are oscillatory. These behaviors are called "transient" and "steady state", respectively. With the symbolic result, you can analyze the result's properties, which is not possible with numeric results.
Visually inspect I1sol and Qsol. They are a sum of terms. Find the terms by using children. Then, find the contributions of the terms by plotting them over [0 15]. The plots show the transient and steady state terms.
I1terms = children(I1sol);
Qterms = children(Qsol);

subplot(1,2,1)
fplot(I1terms,[0 15])
ylim([-2 2])
title('Current terms')

subplot(1,2,2)
fplot(Qterms,[0 15])
ylim([-2 2])
title('Charge terms')
分离稳态和指数项
I1transient = I1terms(has(I1terms,'exp'))
I1steadystate = I1terms(~has(I1terms,'exp'))

Qtransient = Qterms(has(Qterms,'exp'))
Qsteadystate = Qterms(~has(Qterms,'exp'))

%I
subplot(1,2,1)
fplot(I1transient,[0 15])
ylim([-2 2])
title('Current terms')

subplot(1,2,2)
fplot(I1steadystate,[0 15])
ylim([-2 2])
title('Charge terms')

%Q
subplot(1,2,1)
fplot(Qtransient,[0 15])
ylim([-2 2])
title('Current terms')

subplot(1,2,2)
fplot(Qsteadystate,[0 15])
ylim([-2 2])
title('Charge terms')