予選での失態・・・・・・・・・あと固定少数点で微分方程式解けるようになった 2008/10/15(Wed) 01:34:56

動画を取り込んだので、予選での失態を公開




以前ルンゲ・クッタ法での微分方程式を数値解析的に解く方法を試みましたが、計算にかかる時間により一時断念していました。
今回、固定小数点での計算になるよう手直ししてみました。
コレでH8でも制御周期内で次の重心位置の座標を求められるかな？
問題は精度だけどね。

こんな感じのコードです。
1000倍してるだけなんですけどね。

#include <stdio.h>
#include <math.h>

void drunge(long *x, long *y, long *z, long h , long (*f)(long, long, long),long (*g)(long,long,long))
{
 long k1,k2,k3,k4,l1,l2,l3,l4;
 k1 = f(*x, *y, *z);
 l1 = g(*x, *y, *z);
 k2 = f(*x + (h / 2), *y + h * k1 / 2/1000, *z + h * l1 /2/1000);
 l2 = g(*x + h / 2, *y + h * k1 / 2/1000, *z + h * l1 /2/1000);
 k3 = f(*x + h / 2, *y + h * k2 / 2/1000, *z + h * l2 /2/1000);
 l3 = g(*x + h / 2, *y + h * k2 / 2/1000, *z + h * l2 /2/1000);
 k4 = f(*x + h, *y + h * k3/1000, *z + h * l3/1000);
 l4 = g(*x + h, *y + h * k3/1000, *z + h * l3/1000);
 *x += h;
 *y += h * (k1 + 2 * k2 + 2 * k3 + k4)/6/1000;
 *z += h * (l1 + 2 * l2 + 2 * l3 + l4)/6/1000;
 return;
}

long f(long x,long y,long z)
{
 return 0*x + 0*y + z;
}

long g(long x, long y, long z)
{
 return 0*x - 3000*y/1000 - 500*z/1000;
}

int main(void){
 long x = 0, y= 500, z = 1000, xmax = 10000, h;
 int i, n = 100;
 
 h = (xmax - x) / n;
 for(i = 1; i<=n ; i++){
  printf("%4d %d %d %d \n",i,x,y,z);
  drunge(&x,&y,&z,h,f,g);
 }
 return 0;
}

あと、今までルートってどうやって求めればいいんだろう？と思っていたんですが、ちょっと調べて見るとえらく簡単に求められるみたいですね。ちょっと目からウロコ。
これも固定小数点で計算できるようにしました。
これも１００００倍しているだけなんですけどね。（１００００倍だから、返ってくる値は１００倍になります）

int sqrt(long a){
 long b;
 int i;
 
 a = 10000*a;
 b=a;
 for(i=0;i<20;i++){
 b = (b + (a/b))/2;
 }
 return(int(b));
}