プログラマ専用SNS ミクプラ

独学でC言語を勉強しているのですが、連立二次方程式の解き方でつまづいています。。

x^2+y^2=8
(x-1)^2+y^2=5
(x-1)^2+(y-2)^2=1

このような連立二次方程式を解きたいと考えているのですが、普段紙に書いて解いている部分(代入したり約分したり)をプログラムではどのように書いたらいいのかさっぱりなので教えていただけないでしょうか。よろしくお願いします。

二元で式が3個あるのはいいんですか？
さておき、全部円の方程式てすから、中心座標と
円の半径を使えば何とかなりませんか？

コード:

#include <stdio.h>
 
void main(void) {
  
  printf("式1:%(x-a%)^2+%(y-b%)^2²=c\n");
  
  int a,b,c;
  printf("a:");
  scanf("%d",&a);
  
  printf("b:");
  scanf("%d",&b);
  
   printf("c:");
  scanf("%d",&c);
 
 int az,bz,cz;
 
 az=-2*a;
 bz=-2*b;
 cz=a*a+b*b-c;

  printf("式2:%(x-a%)^2+%(y-b%)^2²=c\n");
  
  int a2,b2,c2;
  printf("a:");
  scanf("%d",&a2);
  
  printf("b:");
  scanf("%d",&b2);
  
   printf("c:");
  scanf("%d",&c2);

 int az2,bz2,cz2;
 
 az2=-2*a2;
 bz2=-2*b2;
 cz2=a2*a2+b2*b2-c2;


  printf("式3:%(x-a%)^2+%(y-b%)^2²=c\n");
  
  int a3,b3,c3;
  printf("a:");
  scanf("%d",&a3);
  
  printf("b:");
  scanf("%d",&b3);
  
   printf("c:");
  scanf("%d",&c3);
 
 int az3,bz3,cz3;
 
 az3=-2*a3;
 bz3=-2*b3;
 cz3=a3*a3+b3*b3-c3;
 
 
int d=az-az2;
int e=bz-bz2;
int f=cz-cz2;

int g=az-az3;
int h=bz-bz3;
int i=cz-cz3;

float x=(float)(f*h-i*e)/(g*e-d*h);
float y=(float)(d*i-f*g)/(e*g-d*h);

printf("x=%f\ny=%f",x,y);


}

a,bがない場合は0を打ち込みます

数学を扱うソフトフェアには、主に以下の２つがあります。
　数値処理ソフト：Matlab, Octave
　数式処理ソフト：Mathematica, MAXIMA

数値処理ソフトは、スクリプトに沿って普通の浮動小数点演算に近い計算をするソフト
数式処理ソフトは、1+2*sqrt(2) + 2+3*sqrt(2)を10.07でなく3+5*sqrt(2)と答えてくれるようなソフト

数式処理ソフトの実装は、オブジェクト指向で行えばスマートに行くでしょうが、たぶんコードを書いているうちに
うんざりしてくるのですでにあるソフトを動かして終わりにしたほうがいいです。

---
さて、この問題の場合、自由度がどれだけあるかによってプログラムの難易度（というか行数）が変わってくると
思います。上の３式を一般的に書くと、

コード:

(x-a1)^2 + (y-b1)^2 = c1
(x-a2)^2 + (y-b2)^2 = c2
(x-a3)^2 + (y-b3)^2 = c3

となります。二次方程式のうち、正円だけになり、楕円、放物線や双曲線は除外できるようです。

そして、「超初心者」さんが式の数を気にしていましたが、式が２つでも解なし、式が４つでも解ありの
ケースが考えられます。つまり、２円が交わっていれば解２つ、接していれば１つ、接触していなければ解なしです。

そんでもって、丁寧にそれぞれの円が接触しているかを判別していきます。
で、それはおそらく解析解をもつのでそれ自体を約分するなり数値解を出すなりすればいいと思います。

学生時代に２円の交点座標を計算した記憶がありますが、式の中にsqrt(*)だけでなくarctan(*)もでてきた
ようにおもいました。

私のほうでも２つの円の交点を求めてみました。

コード:

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

typedef struct { double x; double y; } V;

V add(V a, V b) { V c = {a.x+b.x, a.y+b.y}; return c; }
V sub(V a, V b) { V c = {a.x-b.x, a.y-b.y}; return c; }
double inner_product(V a, V b) { return a.x*b.x+a.y*b.y; }
double outer_product(V a, V b) { return a.x*b.y-a.y*b.x; }
V scalar_multiple(double a, V b) { V c = {a*b.x, a*b.y}; return c; }
double norm(V a) { return sqrt(inner_product(a, a)); }


int main() {
        // ２つの円の中心をA, Bとする
        // ２つの円の交点をD, Eとする
        // 直線ABと直線DEの交点をCとする
        V a = {0.0, 0.0};
        V b = {1.0, 2.0};
        double r1 = sqrt(5.0);
        double r2 = 2.0;

        V AB = sub(b, a);
        // 余弦定理より
        double angleCAE = acos( ( r1*r1 + inner_product(AB,AB) - r2*r2 )
                          / ( 2.0 * norm(AB) * r1 ) );
        V tmp1 = {AB.y, -AB.x};
        V normalAB = scalar_multiple(1.0/norm(AB), tmp1);
        V CE = scalar_multiple( r1 * sin(angleCAE), normalAB );
        V AC = scalar_multiple( r1 * cos(angleCAE) / norm(AB), AB );
        V e = add(a, add(AC, CE));
        V d = add(a, sub(AC, CE));
        printf("交点1: (%2f, %2f), 交点2: (%2f, %2f)\n", d.x, d.y, e.x, e.y);

        return 0;
}

三角関数を使わなければ円の交点を求められないかと思ってましたが、
下記URLによればもっと簡単に解けるようです。

http://detail.chiebukuro.yahoo.co.jp/qa ... 1263968897

なので、Reさんのやり方はあってるのかもしれません。確認はしていないです。

２つの円の交点を求めるプログラム　パート２です。
一応厳密解を求めに行っています。おそらくこれが正解ではないかと。

コード:

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

typedef struct { int a; int b; int c; } CIRCLE; // (x-a)^2+(y-b)^2=c
typedef struct { int a; int b; int c; } LINE; // ax+by=c
typedef struct { int a; int b; int c; } QUAD; // ax^2+bx+c=0
typedef struct { int a, b, c, d; } ANSWER; // (a+b*sqrt(c))/d

LINE sub(CIRCLE c1, CIRCLE c2) {
  LINE l;
  l.a = - 2*c1.a + 2*c2.a;
  l.b = - 2*c1.b + 2*c2.b;
  l.c = c1.c - c2.c - c1.a*c1.a - c1.b*c1.b + c2.a*c2.a + c2.b*c2.b;
  return l;
}

QUAD AxPB2(int a, int b) {
  QUAD q;
  q.a = a*a;
  q.b = 2*a*b;
  q.c = b*b;
  return q;
}

// (x-c.a)^2+(y-c.b)^2 = c.c
// (l.a*x-l.a*c.a)^2+(l.a*y-l.a*c.b)^2 = l.a*l.a*c.c
// これにl.a*x=l.c-l.b*yを代入
QUAD ElimX(CIRCLE c, LINE l) {
  QUAD q1 = AxPB2(-l.b, l.c-l.a*c.a);
  QUAD q2 = AxPB2(l.a, -l.a*c.b);
  QUAD q = {q1.a+q2.a, q1.b+q2.b, q1.c+q2.c-l.a*l.a*c.c};
  return q;
}

ANSWER Solve(QUAD q) {
  ANSWER a = {-q.b, 1, q.b*q.b - 4*q.a*q.c, 2*q.a};
  return a;
}

ANSWER SolveX(LINE l, ANSWER y) {
  ANSWER x;
  x.a = y.d*l.c - y.a*l.b;
  x.b = -l.b*y.b;
  x.c = y.c;
  x.d = l.a*y.d;
  return x;
}

int gcd(int a, int b) {
  if (a == 0) return b;
  return gcd(b%a, a);
}

void Print(ANSWER a) {
  int int_part = 1;
  for (int i = (int)sqrt(a.c); i >= 2; i--) {
    if (a.c % (i*i) == 0) {
      int_part = i;
      break;
    }
  }
  a.b *= int_part;
  a.c /= int_part*int_part;
  if (a.d < 0) {
    a.a *= -1;
    a.b *= -1;
    a.d *= -1;
  }
  if (a.c != 1) {
    int gcd1 = gcd(abs(a.a), gcd(abs(a.b), abs(a.d)));
    a.a /= gcd1;
    a.b /= gcd1;
    a.d /= gcd1;
    printf("(%d", a.a);
    if (a.b > 0) printf("+");
    printf("%d*sqrt(%d))/%d", a.b, a.c, a.d);
  } else {
    a.a += a.b;
    int gcd1 = gcd(abs(a.a), abs(a.d));
    a.a /= gcd1;
    a.d /= gcd1;
    if (a.d != 1) printf("%d/%d", a.a, a.d);
    else printf("%d", a.a);
  }
}

int main() {
  CIRCLE c1 = {0, 0, 5}, c2 = {1, 2, 4};
  LINE l = sub(c1, c2);
  ANSWER y1 = Solve(ElimX(c1, l));
  ANSWER y2 = y1;
  y2.b = -y2.b;
  ANSWER x1 = SolveX(l, y1);
  ANSWER x2 = SolveX(l, y2);
  printf("交点1: [ ");
  Print(x1);
  printf(", ");
  Print(y1);
  printf(" ]\n交点2: [ ");
  Print(x2);
  printf(", ");
  Print(y2);
  printf(" ]\n");
  return 0;
}

このような形でプログラムを作ってみたはいいですが、エラーが起こってしまいます。なにか指摘、アドバイスなどありましたら、お願いします。自分でも再び考えてみます。

コード:

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

typedef struct{
	float a;
	float b;
	float c;
}EN;

int main(void){
	EN data[3];

	int i;
	float a[2],b[2],c[2];

	//方程式1〜3に値を代入//
	printf("(x-a)^2+(y-b)^2=cにおいて\n");
	for(i=0;i<3;i++){
		printf("方程式%dのa,b,c の値？\n",i+1);
		scanf("%f%f%f",&data[i].a,&data[i].b,&data[i].c);
	}

	//k*二次方程式+二次方程式=0  k=-1代入して一次方程式(a*x+b*y=c)を求める//

	for(i=0;i<2;i++){
		a[i]=2*data[i].a-2*data[i+1].a;
		b[i]=2*data[i].b-2*data[i+1].b;
		c[i]=(pow(data[i].a,2)+pow(data[i].b,2)-data[i].c)-(pow(data[i+1].a,2)+pow(data[i+1].b,2)-data[i+1].c);
	}

	//二つの直線がでたので、一次連立方程式の解(x,y)を求める(クラメルの公式)//
	float delta,delta_x,delta_y,x,y;
	delta=a[0]*b[1]-b[0]*a[1];
	delta_x=c[0]*b[1]-b[0]*c[1];
	delta_y=a[0]*c[1]-c[0]*a[1];

	x=delta_x/delta;
	y=delta_y/delta;

	printf(" x=%f\n y=%f\n",x,y);

}

akasann さんが書きました:このような形でプログラムを作ってみたはいいですが、エラーが起こってしまいます。なにか指摘、アドバイスなどありましたら、お願いします。

投稿の際のアドバイスとして、「エラーが起こってしまいます。」とだけ書くのではなく、エラーの内容を具体的に書いたほうが親切です。

コンパイルエラーやリンクエラーの場合は、エラーメッセージを(ユーザー名などの個人情報は隠して)そのまま貼る
実行結果がおかしい場合は、入力と期待する動作と観測される動作を書く

といいでしょう。

そうですね。間違ってる部分を示すのも大切ですよね。
なんとなくですが、エラーの理由が分りました。構造体は実体がないというか
int a;
a=a[0].x;
みたいなやり方ができないことに気づいたのでもっと修復してみます。

akasann さんが書きました: 構造体は実体がない

半分あってます。

コード:

typedef struct { float a, b, c; } EN;

ENはintと同じ、「型」であって、変数ではないです。

コード:

  int a;
  EN b;

対して、上記a, bはともに実体をもっています。

ちなみに、クラスの場合はこれをインスタンスと呼ぶことが多いですが、
プリミティブ型(int等)や構造体の場合はあまりそういう
言い方をしないと思います。

akasann さん

上記プログラムの間違いは、プログラム記述の間違いでなく、方式の間違いではないですか？

円の方程式２つからは直線の式が１つしか出ないです。
なぜなら、その直線は求める２つの交点を通っていないといけないからです。

akasann さんが書きました:なんとなくですが、エラーの理由が分りました。構造体は実体がないというか
int a;
a=a[0].x;
みたいなやり方ができないことに気づいたので

確かにint型のaに添字演算子は使えないですが、何が言いたいのでしょうか？
ちなみに、No: 8のコードはコンパイルは通ることを確認しました。

あんどーなつさん、みけCATさん、いろいろアドバイスありがとうございます。
結局のところ、コードの問題だったのかもしれませんが違う部屋？で行ったら実行されました。
出来上がりとしてこのような形になりました。また、アドバイスなどありましたら、ご自由にお願いします。

コード:

#include<stdio.h>
#include<math.h>
 
typedef struct{
    float a;
    float b;
    float c;
}EN;
 
int main(void){
    EN data[3];
 
    int i;
    float a[2],b[2],c[2];
 
    //方程式1〜3に値を代入//
    printf("(x-a)^2+(y-b)^2=cにおいて\n");
    for(i=0;i<3;i++){
        printf("方程式%dのa,b,c の値？\n",i+1);
        scanf("%f%f%f",&data[i].a,&data[i].b,&data[i].c);
    }
 
    //k*二次方程式+二次方程式=0  k=-1代入して一次方程式(a*x+b*y=c)を求める//
 
    for(i=0;i<2;i++){
        a[i]=2*data[i].a-2*data[i+1].a;
        b[i]=2*data[i].b-2*data[i+1].b;
        c[i]=pow(data[i].a,2)+pow(data[i].b,2)-data[i].c-(pow(data[i+1].a,2)+pow(data[i+1].b,2)-data[i+1].c);
    }
 
    //二つの直線がでたので、一次連立方程式の解(x,y)を求める(クラメルの公式)//
	  
	float delta,delta_x,delta_y,x,y;
    delta=a[0]*b[1]-b[0]*a[1];
    delta_x=c[0]*b[1]-b[0]*c[1];
	delta_y=a[0]*c[1]-c[0]*a[1];

	x=delta_x/delta;
	y=delta_y/delta;
 
	printf(" x=%f\n y=%f\n",x,y);
}

akasann さん

ごめんなさい。勘違いしてました。
３つの円の方程式　→　２つの直線の式　→　１つの点
ですね。それなら解けます。

厳密解じゃなくて数値解で解いたのですね。すごいです。
意図を汲むことができなくて申し訳ないです。

コードについてなのですが、短いし、ある程度秩序だっているし、コメントも書いてあるので、
問題ないと思います。

個人的な趣味でいうと、
・EN構造体を定義したならCYOKUSEN構造体もほしい。または構造体を定義しないほうが一貫性がある
・28行目は110文字程あるが、私は80カラムが好み、120カラムがいいという人もいるからこれでもいい

28行目を改行するならば、

コード:

        c[i]=pow(data[i].a,2)+pow(data[i].b,2)-data[i].c
            -(pow(data[i+1].a,2)+pow(data[i+1].b,2)-data[i+1].c);

がいいかと思います。=と-を合わせたり、pと-を合わせたりするまではいいと思います。
マイナス(-)とインデントがこの行が続いていることを教えてくれます。
しかし、人によってはマイナスを行末に置いたりするので、なんで？、と思ったりします。

コード:

        c[i]=pow(data[i].a,2)+pow(data[i].b,2)-data[i].c-
            (pow(data[i+1].a,2)+pow(data[i+1].b,2)-data[i+1].c);

ちなみに k*二次方程式 + 二次方程式 = 0 に k=-1 を代入して、など書いてありますが
元々の問題の出どころとかわかりますか？クラメルの公式を使用されているので
大学生であると推察するのですが...

あんどーなつさん、いろいろアドバイスありがとうございます。私はまだプログラミング初心者ですので、あんどーなつさんのプログラミングも色々と参考にさせてもらっています。説明も分かりやすいですし・・・

k*円の方程式＋円の方程式＝0・・・式(1)
のでどころですが、高校時代に使っていた参考書を参考にしました。

参考書　改正版チャート式　解放と演習　数学Ⅱ+B

上記の参考書のp126当たりを参考にしました。本書によると、式(1)にk=-1を代入することは、x,yの一次式になるようにするためだと考えられます。そうすることによって、二つの円の交点を通る方程式ができあがるそうです。
説明になっているか分りませんが参考になったらうれしいです。

先ほど投稿したコードの問題点としては、円三つが重なってるところ一点なければ実行できないところですが・・

ありがとうございます。
チャート式ですね。あの本は素晴らしいです。
横浜国立大だったか、ax^2 + bx + c = 0を解けという問題があって、
a != 0 と a = 0 のときで場合分けをしなければならないことにショックを受けた記憶があります。

グラフィックス系(GDI, DirectX等)は高校数学の知識が生きるので非常に面白いと思います。
がんばってください。

ごめんごめん、２点だけ。

まず、お褒めをいただきありがとうございます。少し恥ずかしいです。

あと、場合分けなんですよね。そこまでやるとコードが汚れる（見づらくなる）のでやらなかったのですが、
やりたいひとにやらせておけばいいでしょう T(^o^)T

プログラマ専用SNS ミクプラ

C言語での連立二次方程式の解き方について

C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について

Re: C言語での連立二次方程式の解き方について