아두이노와 미분방정식

항상 그렇지만 오늘하는 주제는 재미로 하는 주제 입니다.

인터넷을 보다가 아두이노 328P, 우노보드가 286 컴퓨터보다 빠르다는 기사를 봤습니다.

과연...? 

사실 본인이 컴퓨터 처음공부 할떄는 286,XT 컴퓨터 이전이라....

XT컴퓨터에 도스가 와.... 오락실에서만 하던 테트리스가 흑백이지만. 플로피디스크에 돌아 갈떄 와 했던 기억이....

그래서 한번해 봤습니다.

우노가 미분방정식을 풀수 있는지...?

결론부터 이야기 하면 2차 미분방정식을 푸는데 0.0008,9초 정도 걸리는것보니 정말 빠르네요.

식이 간단한것도 있지만, 그래도 미분방정식이니 재미있게 한번 보죠.

 

두원의 교점을 구하는 문제 입니다. 

사실 제가 아는 방법만해도 많으니, 방법의 소개는 그만 두고, 

제일 쉬운 뉴턴랍슨 방법으로 풀겠습니다.

수치해석을 전문적으로 공부했다면 뉴턴랍슨에 대해서 알수 있지만,

만약 공부하지 않았다면 식의 유도는 약간 어려울수 있습니다.

그래도 이런것도 있다 정도만 알아도 좋죠.

일단 방정식을 유도하고 그방정식을 한번 미분하면 문제는 쉽게 풀립니다.

단 미분이 되지 않는 방벙식은 다른방법이 있습니다.

여기는 수식이 되지 않기 때문에 길게 풀어서 적어 보겠습니다.

 

f(x,y) = (x-a)*(x-a)+(y-b)*(y-b)-c = 0

이게 원의 방정식입니다. 고등학떄 배웠나 중학교떄 배웠나.. ..기물가물...

 

이것을 한번 미분하면 

df(x,y)/dx = 2(x-a)

df(x,y)/dy = 2(y-a)

입니다.

 

총3개의 식만 있으면 문제는 끝입니다.

이제 코딩을해보죠. 

그리고 우리는 문제를 풀기 위해서 

(2,2)의 중심을 가진 지름 4짜리 원과 

(-2,2)의 중심을 가진 지름 4짜리 원 두개를 가지고 문제를 풀어 보죠.

두원의 교점은  0,2가 되겠죠. 답이 겠죠.

아래코딩은 모두 변수 처리 했으니, 값을 변경하면서 풀어 보는것도 재미 있겠습니다.

원의 교점이 두개 일때는 어떻게 합니까 물어보신다면 가르쳐 줄수는 있지만,

재미로하는것이기 떄문에 필요하면 전문적인 수치해석을 공부하라고 권하고 싶습니다.

 

 

#include <math.h>   double ff1(double a, double b, double c, double x, double y) {  return ((x-a)*(x-a) + (y-b)*(y-b) -c) ; }   double df1dx(double a, double b, double c, double x, double y) {  return (2*(x-a) ); }   double df1dy(double a, double b, double c, double x, double y) {  return (2*(y-b) ); }   void newton_raphson() {     int iter =0;     double x,y;     double dt;     double ff[2];     double aa[2][2];     double xx, yy;     double err = 1e-4;     bool conv = false;     char p[35];            x = y = 0;     double a1,b1,c1;     double a2,b2,c2;     a1 = b1=2.0;      a2= -2.0;     b2 = 2.0;     c1=c2=4.0;          while(!conv && iter < 500)     {         ff[0] = ff1(a1,b1, c1, x, y);         ff[1] = ff1(a2,b2, c2, x, y);                  aa[0][0] = df1dx(a1,b1, c1, x, y);         aa[0][1] = df1dy(a1,b1, c1, x, y);         aa[1][0] = df1dx(a2,b2, c2, x, y);         aa[1][1] = df1dy(a2,b2, c2, x, y);                  dt = aa[0][0]*aa[1][1] - aa[0][1]* aa[1][0];                  if( dt != 0.0)         {           xx = ( aa[1][1]* ff[0] - aa[0][1]*ff[1])/dt;           yy = (-aa[1][0]* ff[0] + aa[0][0]*ff[1])/dt;         }         else          {           Serial.println("error");         }         x -= 0.5*xx;         y -= 0.5*yy;         //sprintf(p, "iter = %d, %s, %s", iter, String(xx,5).c_str(), String(err ,5).c_str());            //Serial.println(p);                 conv = (abs(xx) <= err*x ) && (abs(yy) <= err*y);         iter++;     }       sprintf(p, "iter = %d, %s, %s", iter, String(x,2).c_str(), String(y,2).c_str());        Serial.println(p);      return 0; }   void setup()  {   Serial.begin(9600);   while (!Serial)    {     ; // wait for serial port to connect. Needed for native USB port only   }   Serial.println("initialization done.");   }   void loop()  {     unsigned long currentMillis = millis();          newton_raphson();       Serial.print("duration time= ");     Serial.println(millis() -currentMillis);     Serial.println();     Serial.println();     delay(5000); }   Colored by Color Scripter