이 번역 페이지는 최신 내용을 담고 있지 않습니다. 최신 내용을 영문으로 보려면 여기를 클릭하십시오.
라플라스 변환을 사용하여 미분 방정식 풀기
Symbolic Math Toolbox™에서 다음 워크플로로 라플라스 변환을 사용하여 미분 방정식을 풉니다. 라플라스 변환의 간단한 예제는 laplace 및 ilaplace 항목을 참조하십시오.
정의: 라플라스 변환
함수 f(t)의 라플라스 변환은 다음과 같습니다.
F (s)=∫0∞f(t)e-tsdt.
개념: 기호 워크플로 사용하기
기호 워크플로에서는 계산을 수치 형식이 아닌 자연적인 기호 형식으로 유지합니다. 이 접근 방식은 해의 속성을 이해하고 정확한 기호 값을 사용하게 해줍니다. 수치 결과가 필요할 때나 기호적으로 진행할 수 없을 때만 기호 변수에 숫자를 대입합니다. 자세한 내용은 수치 연산방식 또는 기호 연산방식 선택하기 항목을 참조하십시오. 일반적으로 단계는 다음과 같습니다.
방정식을 선언합니다.
방정식을 풉니다.
-
값을 대입합니다.
결과를 플로팅합니다.
결과를 분석합니다.
워크플로: 라플라스 변환을 사용하여 RLC 회로 풀기
방정식 선언하기
라플라스 변환을 사용하여 초기 조건이 있는 미분 방정식을 풀 수 있습니다. 예를 들어, 다음 회로와 같은 저항-인덕터-커패시터(RLC) 회로를 풀 수 있습니다.
저항(단위: 옴): R1,R2,R 3
전류(단위: 암페어): I1,I2,I3
유도용량(단위: 헨리): L
정전용량(단위: 패럿): C
AC 전압원(단위: 볼트): E(t)
커패시터 전하(단위: 쿨롬): Q(t )
키르히호프의 전압 및 전류 법칙을 적용하여 다음 방정식을 얻습니다.
I1=I2+ I3LdI1dt+I1R1+I2R2=0E(t)+I2R2- QC-I3R3=0
위의 방정식에 관계식 I3=dQ/dt(커패시터 충전 속도)를 대입하여 RLC 회로의 미분 방정식을 얻습니다.
dI1 dt-R2LdQdt=-R1+R2LI1dQdt=1R2+R3(E(t)-QC)+R2R 2+R3I1
변수를 선언합니다. 물리량이 양수 값을 가지므로, 변수에 이에 해당하는 가정을 설정합니다. E (t)를 1V의 교류 전압이라고 하겠습니다.
syms L C I1(t) Q(t) s R = sym('R%d',[1 3]); assume([t L C R] > 0) E(t) = 1*sin(t); % AC voltage = 1 V
미분 방정식을 선언합니다.
dI1 = diff(I1,t); dQ = diff(Q,t); eqn1 = dI1 - (R(2)/L)*dQ == -(R(1)+R(2))/L*I1
eqn1(t) =
∂∂t I1(t)-R2 ∂∂t Q(t)L=-I1(t) R1+R2L
eqn2 = dQ == (1/(R(2)+R(3))*(E-Q/C)) + R(2)/(R(2)+R(3))*I1
eqn2(t) =
∂∂t Q(t)=sin(t)-Q(t)CR2+R3+R2 I1(t)R2+R3
방정식 풀기
eqn1과 eqn2의 라플라스 변환을 계산합니다.
eqn1LT = laplace(eqn1,t,s)
eqn1LT =
s laplace(I1(t),t,s)-I1(0)+R2 Q(0)-s laplace(Q(t),t,s)L=-R1+R2 laplace(I1(t),t,s)L
eqn2LT = laplace(eqn2,t,s)
eqn2LT =
s laplace(Q(t),t,s)-Q(0)=R2 laplace(I1(t),t,s)R2+R3+Cs2+1-laplace(Q(t),t,s)C R2+R3
함수 solve는 기호 변수에 대해서만 방정식을 풉니다. 따라서 solve를 사용하려면 먼저 laplace(I1(t),t,s)와 laplace(Q(t),t,s)에 변수 I1_LT와 Q_LT를 대입하십시오.
syms I1_LT Q_LT eqn1LT = subs(eqn1LT,[laplace(I1,t,s) laplace(Q,t,s)],[I1_LT Q_LT])
eqn1LT =
I1,LT s-I1(0)+R2 Q(0)-QLT sL=-I1,LT R1+R2L
eqn2LT = subs(eqn2LT,[laplace(I1,t,s) laplace(Q,t,s)],[I1_LT Q_LT])
eqn2LT =
QLT s-Q(0)=I1,LT R2R2+R3-QLT-Cs2+1C R2+R3
I1_LT와 Q_LT에 대해 방정식을 풉니다.
eqns = [eqn1LT eqn2LT]; vars = [I1_LT Q_LT]; [I1_LT, Q_LT] = solve(eqns,vars)
I1_LT =
L I1(0)-R2 Q(0)+C R2 s+L s2 I1(0)-R2 s2 Q(0)+C L R2 s3 I1(0)+C L R3 s3 I1(0)+C L R2 s I1(0)+C L R3 s I1(0)s2+1 R1+R2+L s+C L R2 s2+C L R3 s2+C R1 R2 s+C R1 R3 s+C R2 R3 s
Q_LT =
C R1+R2+L s+L R2 I1(0)+R1 R2 Q(0)+R1 R3 Q(0)+R2 R3 Q(0)+L R2 s2 I1(0)+L R2 s3 Q(0)+L R3 s3 Q(0)+R1 R2 s2 Q(0)+R1 R3 s2 Q(0)+R2 R3 s2 Q(0)+L R2 s Q(0)+L R3 s Q(0)s2+1 R1+R2+L s+C L R2 s2+C L R3 s2+C R1 R2 s+C R1 R3 s+C R2 R3 s
I1_LT와 Q_LT의 라플라스 역변환을 계산하여 I1과 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/4F, L=1 .6H라고 하겠습니다. 초기 전류는 I1(0)= 2A이고 초기 전하는 Q(0)=2C라고 가정합니다.
vars = [R L C I1(0) Q(0)]; values = [4 2 3 1.6 1/4 2 2]; I1sol = subs(I1sol,vars,values)
I1sol =
200 cos(t)8161+405 sin(t)8161+16122 e-81 t40 cosh(1761 t40)-742529 1761 sinh(1761 t40)141954218161
Qsol = subs(Qsol,vars,values)
Qsol =
924 sin(t)8161-1055 cos(t)8161+17377 e-81 t40 cosh(1761 t40)+1109425 1761 sinh(1761 t40)306008978161
결과 플로팅하기
전류 I1sol과 전하 Qsol을 플로팅합니다. 2개의 서로 다른 시간 구간 0≤t≤15와 2≤t≤25를 사용하여 과도 상태 동작과 정상 상태 동작을 모두 표시합니다.
subplot(2,2,1) fplot(I1sol,[0 15]) title('Current') ylabel('I1(t)') xlabel('t') subplot(2,2,2) fplot(Qsol,[0 15]) title('Charge') ylabel('Q(t)') xlabel('t') subplot(2,2,3) fplot(I1sol,[2 25]) title('Current') ylabel('I1(t)') xlabel('t') text(3,-0.1,'Transient') text(15,-0.07,'Steady State') subplot(2,2,4) fplot(Qsol,[2 25]) title('Charge') ylabel('Q(t)') xlabel('t') text(3,0.35,'Transient') text(15,0.22,'Steady State')
결과 분석하기
초기에는 전류와 전하가 지수적으로 감소합니다. 그러나 장기적으로는 진동합니다. 이러한 동작은 각각 "과도 상태"와 "정상 상태"라고 합니다. 기호 결과를 사용하면 결과의 속성을 분석할 수 있습니다. 이는 수치 결과를 사용할 때는 가능하지 않습니다.
I1sol과 Qsol을 시각적으로 검토합니다. 항들의 합을 볼 수 있습니다. children을 사용하여 항을 구합니다. 그런 다음 [0 15]에 대해 항을 플로팅하여 항의 비중을 구합니다. 플롯에서 과도 상태 항과 정상 상태 항을 볼 수 있습니다.
I1terms = children(I1sol); I1terms = [I1terms{:}]; Qterms = children(Qsol); Qterms = [Qterms{:}]; figure; subplot(1,2,1) fplot(I1terms,[0 15]) ylim([-0.5 2.5]) title('Current terms') subplot(1,2,2) fplot(Qterms,[0 15]) ylim([-0.5 2.5]) title('Charge terms')
플롯을 통해 I1sol에는 하나의 과도 상태 항과 하나의 정상 상태 항이 있고 Qsol에는 하나의 과도 상태 항과 두 개의 정상 상태 항이 있음을 볼 수 있습니다. 시각적 검토를 통해 I1sol과 Qsol에 exp 함수를 포함하는 하나의 항이 있는 것을 알 수 있습니다. 이 항이 과도 지수 감쇠의 원인이 된다고 가정합니다. has를 사용하여 exp가 있는 항이 있는지 확인함으로써 I1sol과 Qsol의 과도 상태 항과 정상 상태 항을 분리합니다.
I1transient = I1terms(has(I1terms,'exp'))
I1transient =
16122 e-81 t40 cosh(1761 t40)-742529 1761 sinh(1761 t40)141954218161
I1steadystate = I1terms(~has(I1terms,'exp'))
I1steadystate =
(200 cos(t)8161405 sin(t)8161)
마찬가지로, Qsol을 과도 상태 항과 정상 상태 항으로 분리합니다. 이 결과는 기호 계산이 문제를 분석하는 데 어떻게 도움이 되는지 보여줍니다.
Qtransient = Qterms(has(Qterms,'exp'))
Qtransient =
17377 e-81 t40 cosh(1761 t40)+1109425 1761 sinh(1761 t40)306008978161
Qsteadystate = Qterms(~has(Qterms,'exp'))
Qsteadystate =
(-1055 cos(t)8161924 sin(t)8161)