매트랩 [2계 상미분방정식]의 연립 방정식으로의 변환 (수치적, 해석적 풀이)

수치적 풀이와 해석적 풀이를 쉽게 설명하면, 수치적 풀이는 미분을 하듯, 실수 구간을 쪼개서 각각의 값을 대입하여 푸는 방법이고, 해석적 풀이는 어떠한 식의 공식을 구하는 방법이다. 

 

2계 상미분방정식을 그냥 풀 수도 있지만, 연립 방정식으로 변환 후 풀 수 있다. 이때의 장점은 2계가 1계로 변한다는 것이다.

 

감쇠운동을 하는 용수철의 진동에 대해 풀어보자.

감쇠운동하는 용수철의 진동 공식 :             질량: m ,  감쇠상수 : c ,  용수철 상수 : k , 용수철 변위 : y

                                                                              m * y'' = - c * y' - k * y

 

 

1. 먼저 상수를 설정해주자.

% mass of spring
m = 1;

% damping constant
c = 2;

% spring constant
k = 0.75;

2. 용수철 공식을 matlab 안에서 함수로 나타내주자. 이때 matlab 안에 n계 상미분 방정식을 n 개의 연립상미분방정식으로 바꿔주는 odeToVectorField를 사용하자.

syms y(t)
[V] = odeToVectorField(m*diff(y,2) == -c*diff(y)-k*y);

V의 class 는 'sym' 이다. 즉 우리가 손으로 쓰는 일반적인 수학 식이다.

Y[1] = y 이고 Y[2] = y' 이다. 

V =                 Y[2]                         :  첫번째 식은 Y[1]' = Y[2]                         
                      - (3*Y[1])/4 - 2*Y[2]  : 두번째 식은 Y[2]' = - (3*Y[1])/4 - 2*Y[2] (즉, 감쇠운동하는 용수철의 진동 공식이다)

 

 

3. 수치해석으로 먼저 풀어보자.

수치해석을 위한 matlab의 함수 중 ode45 를 사용하려면 class가 sym 이 아닌 function 이 와야 한다. 따라서 V 의 class를 function handle 을 {'t', 'Y'} 로 갖고 있는 function으로 바꿔주자. 그리고 ode45 를 통해 풀어준다.

M = matlabFunction(V, 'vars', {'t', 'Y'});

% [0 20] : x축 범위 , [1;1] : 초기값 설정 Y[1](0) = 1; Y[2](0) = 1
[X, Y] = ode45(M, [0 20], [1; 1]);

X 안에는  0-20까지 수가 총 81개 들어 있다. (81x1 matrix)

X =

         0
    0.0183
    0.0365
    0.0548
    0.0731
    0.1644
    0.2558
    0.3471
    0.4384
    0.5496
    0.6609
    0.7721
    0.8833
    0.9768
    1.0702
    1.1637
    1.2572
    1.3753
    1.4933
    1.6114
    1.7295
    1.8627
    1.9960
    2.1293
    2.2625
    2.4152
    2.5679
    2.7205
    2.8732
    3.0479
    3.2227
    3.3974
    3.5722
    3.7723
    3.9725
    4.1726
    4.3727
    4.6044
    4.8361
    5.0677
    5.2994
    5.5733
    5.8473
    6.1212
    6.3951
    6.7345
    7.0740
    7.4134
    7.7528
    8.1934
    8.6340
    9.0745
    9.5151
    9.9191
   10.3231
   10.7270
   11.1310
   11.5327
   11.9344
   12.3361
   12.7378
   13.1348
   13.5318
   13.9289
   14.3259
   14.7209
   15.1159
   15.5109
   15.9059
   16.2999
   16.6940
   17.0880
   17.4820
   17.8857
   18.2894
   18.6932
   19.0969
   19.3227
   19.5484
   19.7742
   20.0000

Y에는 X의 각 숫자를 Y[1] (y) 와 Y[2] (y')에 대입안 값이 각각 1열, 2열에 들어있다.

Y =

    1.0000    1.0000
    1.0178    0.9505
    1.0347    0.9026
    1.0508    0.8562
    1.0660    0.8113
    1.1306    0.6068
    1.1779    0.4331
    1.2105    0.2859
    1.2308    0.1617
    1.2416    0.0368
    1.2400   -0.0634
    1.2283   -0.1430
    1.2087   -0.2056
    1.1875   -0.2471
    1.1628   -0.2801
    1.1353   -0.3058
    1.1058   -0.3253
    1.0663   -0.3425
    1.0252   -0.3529
    0.9832   -0.3578
    0.9408   -0.3583
    0.8933   -0.3549
    0.8464   -0.3481
    0.8006   -0.3388
    0.7562   -0.3277
    0.7072   -0.3136
    0.6605   -0.2984
    0.6161   -0.2827
    0.5742   -0.2669
    0.5291   -0.2490
    0.4871   -0.2316
    0.4481   -0.2149
    0.4120   -0.1989
    0.3739   -0.1817
    0.3392   -0.1657
    0.3075   -0.1509
    0.2787   -0.1372
    0.2486   -0.1228
    0.2217   -0.1098
    0.1976   -0.0981
    0.1761   -0.0875
    0.1537   -0.0765
    0.1341   -0.0668
    0.1170   -0.0583
    0.1020   -0.0509
    0.0861   -0.0430
    0.0727   -0.0363
    0.0614   -0.0307
    0.0518   -0.0259
    0.0415   -0.0208
    0.0333   -0.0167
    0.0267   -0.0134
    0.0215   -0.0107
    0.0175   -0.0088
    0.0143   -0.0072
    0.0117   -0.0059
    0.0096   -0.0048
    0.0078   -0.0039
    0.0064   -0.0032
    0.0052   -0.0026
    0.0043   -0.0021
    0.0035   -0.0018
    0.0029   -0.0014
    0.0024   -0.0012
    0.0019   -0.0010
    0.0016   -0.0008
    0.0013   -0.0007
    0.0011   -0.0005
    0.0009   -0.0004
    0.0007   -0.0004
    0.0006   -0.0003
    0.0005   -0.0002
    0.0004   -0.0002
    0.0003   -0.0002
    0.0003   -0.0001
    0.0002   -0.0001
    0.0002   -0.0001
    0.0002   -0.0001
    0.0001   -0.0001
    0.0001   -0.0001
    0.0001   -0.0001

4. 그래프 그리기.

x축의 값은 X안에 들어 있고, 우리가 y축의 값으로 필요한 것은 Y의 1열(Y[1])이다.

plot(X,Y(:,1), 'o');

수치해석으로 그린 감쇠-용수철 진동

 

5. 검증을 위해 해석적으로 풀어보자.

syms g(x)
ode = m*diff(g,x,2) == -c*diff(g) -k*g;
Dg = diff(g);

% 초기값 조건 설정
cond1 = g(0) == 1;
cond2 = Dg(0) == 1;

% 2계 상미분방정식 풀기
ySol = dsolve(ode, [cond1 cond2]);

6. 그래프 그리기.

subplot(2,1,2);
hold on
fplot(ySol, [0 20]);
hold off

이상으로 검증을 마친다!