내가 공부한 옥타브 22.방향장

매트랩을 이용하여 제1계 상미분 방정식의 방향장을 그려보는 연습을 해보겠습니다.

방향장? 그 각점 x,y에서 미분방정식의 해가 있을텐데 그 그래프의 해의 접선을 다 구했을때 그 접선의 기울기를 각 점마다 나타낸것을기울기의 장이다. 다른 말로는  방향장이라고 합니다.

방향장을 보면 해의 모양을 이은것 같은 모양이라는 것을 추측할 수 있습니다. 

 

방향장을 x,y평면에 어떻게 나타내냐

접선을 그리되 겹치지 않도록 아주 짧은 부분만 (선소)그려서 쭉 나열하면 그것이 방향장이 됩니다. 

y

전기장?

장?

 Field

 

이런 아이디어를 구현하려면 손으로 일일히 그리기가 힘들기 때문에, 프로그램을 이용하여 그리는것이 오늘의 목표입니다.

 

 

 

매트랩함수 중에서quiver이라는 함수가 있는데 이 함수를 이용해서 방향장을 그려보고자 합니다. 

사용방법을 익혀보도록 하겠습니다. 

아래 그림에서 보다시피 4개의 인자를 가지고 있는데, X,Y는 XY평면에 각 점을 찍어서 점의 좌표를 제공하고 

U,V는 그 점에 그릴 벡터의 X,Y성분을 제공하면 그 지점의 방향장에 해당되는 각 점의 벡터를 그리게 됩니다. 

일단 X,Y평면의 좌표를 줍니다. 

대문자 X,대문자Y에 좌표를 집어넣습니다. 

좌표를 집어넣는 방법은 행렬로 대입해야하는데 하나하나 다 찍어낼 수 없으니까,

10과에서 3차원 그래프를 공부할때 meshgrid() 라는 함수를 공부했습니다.

meshgrid()함수는 X의 범위를 정해주고  Y의 범위를 정해주면 그것에 따라서 X,Y좌표로 꽉 채워넣게 됩니다.

 

 

([X의 범위]  [Y의 범위])

시작값: 증분: 끝값

X,Y각점에 대해서 X성분 U, Y성분 V를 주어야 하는데 , U는 기본적으로 1이라는 값을 주면 V값만 상대적인 값을 주면

벡터로 그 상대적인 크기가 나온다.

 

U에 전부 1을 채우려면, ones()라는 함수를 사용하면 됩니다. 

괄호안에는 행렬의  size를 입력하면 됩니다. 

이렇게 하면 좌표 X,Y 그 지점에 그릴 벡터 U,V 의 값을 다 주었기 때문에 quiver()함수를 사용해서 

다음과 같이 나타낼 수 있습니다. 

 

전체 코드입니다. 

[X Y]=meshgrid( [-5:0.5:5],[-5:0.5:5] );
U=ones(size(X));
V=U; 
quiver(X,Y,U,X) %X,Y좌표가 지정하는 점에 벡터 U,V를 그려라.

 

 

실제로 겹치지는 않지만 벡터의 길이가 너무 길어서 마음에 안든다면 조금 변경시킬 수도 있습니다. 

인자를 하나더 넣어서 벡터의 크기를 조절할 수 있습니다. 

이러한 벡터를 scaling vector이라고 합니다. 

전체코드입니다. 

[X Y]=meshgrid( [-5:0.5:5],[-5:0.5:5] );
U=ones(size(X));
V=U; 
quiver(X,Y,U,V,0.5) %X,Y좌표가 지정하는 점에 벡터 U,V를 그려라.

 

 

 

똑같이 X,Y평면의 각 점에 벡터를 그리는데 X가1보다 클때는 벡터 1,2를 그리지만

 

X가 -2와 1사이에서는

 

역시 똑같이 하면 됩니다 .

그런데,

 

X좌표와 Y좌표를 찾아내는데 meshgrid로 만들어내겠다. 

 

 

벡터[1,1]을 그리자고 했으니, V도 역시 모두1로 채우면 됩니다. 

좌표  X,Y랑 그 좌표에 그릴 벡터 U,V를 다 주었기 때문에  quiver()이라는 함수를 사용해서 x,y좌표가 지정하느 ㄴ점에 벡터UV를 그리라고 하면 다음과 같이 나타납니다. 

 

전체코드입니다. 

%18a1-2번 입니다. 

[X Y]=meshgrid([-5:5],[-5:5]);
U=ones(size(X)); %U를 X의 사이즈 만큼 모두 1로 채운다. 
V=U
quiver(X,Y,U,V);

 

만약 벡터의 길이가 너무길어서 마음에 안든다면, 또는 길어서 서로 겹친다면 scaleling vector을 추가로 작성하면 됩니다. 

전체코드와 Figure창입니다.

%18a1-2번 입니다. 

[X Y]=meshgrid([-5:5],[-5:5]);
U=ones(size(X)); %U를 X의 사이즈 만큼 모두 1로 채운다. 
V=U;
quiver(X,Y,U,V,0.5);

 

 

 

간격을 늘려서 그리는 점의 개수가 조금 줄어들었습니다.

 

 

 

 

U에다가는 X의 사이즈 만큼 집어넣고, V에다가는 값을 넣어야 되는데 모두 1을 넣는게 아니라 X의 범위에 따라 -1 0 1을 넣어야 하므로

※범위에 따라 서로다른 값을집어넣는 방법

 

전체코드와 결과화면 입니다 

 

%ex18a1_3
[X Y]=meshgrid([-5:5],[-5:5]); %간격을 1로
U=ones(size(X));
V=(X<-1).*(-1)+(X>1).*(1)+(X>-1&X<1).*0    %마지막에 0을 곱한항은 작성하지 않아도 됩니다. 
quiver(X,Y,U,V)

 

 

 

quiver()함수를 사용할 줄 알았으니 이 함수를 사용해서 주어진 미분 방정식의 방향장을 그려보도록 하겠습니다. 

 

전체코드입니다. 

% ex18a2_1
[X,Y]=meshgrid([-5:0.5:5], [-5:0.5:5]);
U=ones(size(X));
V=0.2.*X.*Y;
quiver(X,Y,U,V)

 

위그림에서 수정하고 싶은부분: 벡터의크기가 일정하지않다. 

크기를 일정하게 하고 싶으면 단위벡터 (크기가 1인 벡터)를 집어넣으면 되겠습니다. 

 

 

 

각 벡터의 크기를 구합니다. 

(수학 공식)

다음과 같은 코드를 추가하면 벡터의 크기를 구할 수 있게됩니다. 

L=sqrt(X.^2+Y.^2)

 

자기자신의 크기로 다 나누어주면  크기가 거의 같은 다음의 Filed를 완성하게 됩니다. 

% ex18a2_2
[X,Y]=meshgrid([-5:0.5:5], [-5:0.5:5]);
U=ones(size(X));
V=0.2.*X.*Y;
L=sqrt(X.^2+Y.^2)
quiver(X,Y,U./L,V./L)

 

중심부분이 저렇게 나타난 이유는 XY가 0의 근처라서 계산의 오류가 나서 저렇게 된것입니다.

 

바깥의 여백들이 마음에 안들어서 axis tight() 함수를 사용해서 

여백을 채워주면 됩니다. 

 

결과화면 입니다. 

 

전체 코드입니다. 

% ex18a2_3
[X,Y]=meshgrid([-5:0.5:5], [-5:0.5:5]);
U=ones(size(X));
V=0.2.*X.*Y;
L=sqrt(X.^2+Y.^2)
axis tight
quiver(X,Y,U./L,V./L,'r')

 

벡터장의 한 가운데 (0,0)좌표에서 벡터의 크기가 조금 다른데 저 부분을 손보는 방법이 있습니다. 

 

이제 우리가 예제 18b에서는 다른 1계선형미분방정식의 방향장을 다른 범위에 대해서 그려봅니다. 

위에서 했던 것 과 같이 같은 흐름의 코드로 작성하면 됩니다. 

 

방향장의 모양이 해의 모양을 대체로 보여준다고 생각하면 됩니다. 

 

 

그 다음에는 방향장 뿐만아니라 contour flot도 같이그리는 연습을 해보겠습니다. 

contour flot 이라는 것은 어떤 함수가 있을떄 z=f(x,y)이렇게 변수 2개 짜리함수가 있을때 

z값이 일정한 점들을 이어서 만든 plot입니다. 지도에서 봤을떄 높이가 똑같은 점들끼리 연결해서 만든 등고선, 

함수에서 봤을때는 함수값이 똑같은 선만 이어서 만든 등치선, 이와 같은 것들을 contour flot라고 부릅니다. 

 

contour를 그리는 방법은 X,Y평면에 점으 ㄹ찍어야 하므로, X,Y좌표를 먼저 찍어주면 됩니다. 

좌표를 정했으면 그 좌표에 대한 함숫값을 계산해주어야 합니다. 

전체코드입니다. 

%%ex18c1_1

[X Y]=meshgrid([-5:0.2:5],[-2:0.2:2]);
Z=X.^2+Y.^4;
contour(X,Y,Z)

 

 

자동으로 C값을 바꾸어 가면서 그리는 모양입니다. 

 

 

 

만약 C값을 정해서 그리고 싶다고 하면 contour의 마지막 인자에 C값의 벡터를 설정해주면 됩니다. 

%%ex18c1_1

[X Y]=meshgrid([-5:0.2:5],[-2:0.2:2]);
Z=X.^2+Y.^4;
contour(X,Y,Z,[1,4,16])

 

 

Figure 화면 입니다. 

 

 

방향장이 해를 보여준다는 것을 알고 있으므로 다음예제에서 주어진 1계미분방정식에서 방향장도 그리고 , 그 위에 겹쳐서 contour도 같이 그려보라는 예제입니다. 

방향장을 그리는데에는 주어진 미분방정식만 있으면 되지만 contour를 그리기 위해서는 우리가 이것의 해를 알아야 합니다.

해를 직접 분리형 미분방정식을 풀어서 구했습니다.

이해를 직접 등치선으로 그려보도록 하겠습니다.

 

 

meshgrid로 좌표를 지정해주고

 

 

 

전체코드

 

왜 등치선이 나오는 거지?

 

 

 

전체코드 와 Figure창 입니다. 

%ex18c2_1
[X Y]=meshgrid([-5:0.2:5],[-2:0.2:2]);
Z=X.^2+Y.^4;
contour(X,Y,Z);
contour(X,Y,Z,[1 4 16]);

[X Y]=meshgrid([0:0.5:15],[-4:0.5:4]);
U=ones(size(X));
V=(4-X)./(Y.^3+2) %기울기
L=sqrt(U.^2+V.^2);
quiver(X,Y,U./L,V./L,0.5,'r')
axis tight

Z=Y.^4+8*Y-16*X+2*X.^2
hold on
contour(X,Y,Z)

 

 

 

 

방향장을 이은 곡선이 해의 그래프를 나타낸다는 것을 알아보았습니다.