라플라스 미분방정식 계산기 - lapeullaseu mibunbangjeongsig gyesangi

라플라스 변환을 사용하여 미분 방정식 풀기

Symbolic Math Toolbox™에서 다음 워크플로로 라플라스 변환을 사용하여 미분 방정식을 풉니다. 라플라스 변환의 간단한 예제는 laplace 및 ilaplace 항목을 참조하십시오.

정의: 라플라스 변환

함수 f(t)의 라플라스 변환은 다음과 같습니다.

F (s)=∫0∞f(t)e-tsdt.

개념: 기호 워크플로 사용하기

기호 워크플로에서는 계산을 수치 형식이 아닌 자연적인 기호 형식으로 유지합니다. 이 접근 방식은 해의 속성을 이해하고 정확한 기호 값을 사용하게 해줍니다. 수치 결과가 필요할 때나 기호적으로 진행할 수 없을 때만 기호 변수에 숫자를 대입합니다. 자세한 내용은 수치 연산방식 또는 기호 연산방식 선택하기 항목을 참조하십시오. 일반적으로 단계는 다음과 같습니다.

1. 방정식을 선언합니다.

2. 방정식을 풉니다.

3. 값을 대입합니다.

4. 결과를 플로팅합니다.

5. 결과를 분석합니다.

워크플로: 라플라스 변환을 사용하여 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)