매트랩 [1계 상미분방정식] 초기값 문제. 기울기 벡터장 그리기.

공학적인 문제를 풀 때의 과정은

물리 시스템 - 수학적 모델 - 수학적 풀이 - 물리적 해석을 거친다. 이때 우리가 컴퓨터의 도움을 효율적으로 받을 수 있는 부분이 수학적 풀이 부분이다. 모델까지는 다양한 변수가 포함될 수 있기에 사람의 판단력이 필요하다고 생각한다.

 

간단한 1계 상미분 방정식 초기값 문제 y' = x+y, y(0) = 1 이 주어졌을 때 우리는 좌표 (x, y)에서 기울기 y' = x + y

를 갖는 함수를 떠올릴 수 있다. 그렇다면 함수의 기울기장을 그리면 함수의 그래프의 개형을 알 수 있지 않을까?

 

전체 코딩

syms y(x)
ode = diff(y,x) == x+y;
cond = y(0) == 1;
ySol(x) = dsolve(ode, cond);

[a , b] = meshgrid(-2:0.2:3, -1:0.2:2);
s = a + b;
L = sqrt(1 + s.^2);

quiver(a,b, 1./L, s./L, 0.5), axis tight
xlabel 'x', ylabel 'y'
title 'Direction field for dy/dx = x+y'

hold on
fplot(ySol, [-2,3])
hold off

 

1. 상미분 방정식 설정. 주어진 문제에 알맞는 상미분 방정식을 정의. (ode = ordinary differential equation. 상미분 방정식)

syms y(x)
ode = diff(y,x) == x+y;

2. 기울기장(벡터장)을 표시할 좌표 공간 설정. meshgrid를 통해 x와 y에는 좌표공간에서의 x좌표와 y좌표가 각각 존재.

    기울기 벡터장 s를 정의해 준다.

[a , b] = meshgrid(-2:0.2:3, -1:0.2:2);
s = a + b;

 3. 기울기장 그리기. quiver를 이용. axis = tight 은 그래프에 기울기장이 딱 알맞게 들어간다.

quiver(a, b, ones(size(s)),s), axis tight

4. 기울기장의 크기를 일정하게 하기 위해 모든 기울기장의 크기를 1로 통일하자.

L = sqrt(1 + s.^2);
quiver(a,b, 1./L, s./L, 0.5), axis tight
xlabel 'x', ylabel 'y'
title 'Direction field for dy/dx = x+y'

결과:

자 이제 이 기울기장을 검증해 보자!

아직 안 쓴 초기 조건 y(0) = 1을 사용하여 상계미분방정식을 풀어서 ySol에 저장해 주자.

syms y(x)
ode = diff(y,x) == x+y;
cond = y(0) == 1;
ySol(x) = dsolve(ode, cond);

ySol =

이 나오고 이를 벡터장과 같은 그래프에 그려보자.

hold on/off 를 통해 같은 figure 위에 두 그래프가 공존하게 된다.

L = sqrt(1 + s.^2);
quiver(a,b, 1./L, s./L, 0.5), axis tight
xlabel 'x', ylabel 'y'
title 'Direction field for dy/dx = x+y'
hold on
fplot(ySol, [-2,3])
hold off

따라서 초기조건을 적용하여 상미분 방정식을 풀었을 때 기울기장을 따라가는 그래프가 그려짐을 확인할 수 있다!