掲示板利用宣言

 次のフォームをすべてチェックしてからご利用ください。

 私は

 題名と投稿者名は具体的に書きます。
 課題の丸投げはしません。
 ソースの添付は「HTML変換ツール」で字下げします。
 返信の引用は最小限にします。
 環境(OSとコンパイラ)や症状は具体的に詳しく書きます。
 返信の付いた投稿は削除しません。
 マルチポスト(多重投稿)はしません。

掲示板2

管理者用メニュー    ツリーに戻る    携帯用URL    ホームページ    ログ    タグ一覧

No.29218

プログラムニュートンラフソン法
投稿者---ヒロ(2006/12/12 23:34:18)


下記の非線型連立方程式を

f1(x1,x2)=x1^2+x2^2-1=0
f2(x1,x2)=(x1+x2)^2-4x1-1=0

ニュートンラフソン法を使ってときたいのですが、
ヤコビ行列やらなんやらでわけがわからなくなりました。
一応ですが、自分なりにプログラムを作りましたが未完成です。
どうすればいいでしょうか?ご教授おねがいします!!

#include<iostream>
#include<cmath>
#include<iomanip>
using namespace std;

//前準備

double f(double x1,double x2)  {return x1*x1+x2*x2-1;}

double dfx(double x1,double x2){return 2*x1;}

double dfy(double x1,double x2){return 2*x2;}

double g(double x1,double x2)  {return pow(x1+x2,2)-4*x1-1;}

double dgx(double x1,double x2){return (2*x1+2*x2-4);}

double dgy(double x1,double x2){return 2*x1+2*x2;}

double J(double x1,double x2){    //逆行列
  double k;
  k=1/(4*(x1*x1-x2*x2+2*x2));
  return (dfx(x1,x2)*dgy(x1,x2)-dfy(x1,x2)*dgx(x1,x2));
}

main(){
  int k=0;
  double x,y,x1,y1;
  double eps=1e-6;
  x1=1.0;
  y1=1.0;

   do{
     x=x1-(f(x1,y1)*dgy(x1,y1)-g(x1,y1)*dfy(x1,y1))/J(x1,y1);
     y=y1-(g(x1,y1)*dfx(x1,y1)-f(x1,y1)*dgx(x1,y1))/J(x1,y1);

     x1=x;
     y1=y;
     k++;

     cout<<showpoint<<"x="<<setw(13)<<x1<<setw(5)<<" y= "<<y1<<endl;
	 
      if(k>100){
	cout<<"error!!!"<<endl;
         return 1;
      }

   }while(fabs(f(x1,y1)-g(x1,y1))>eps);


  return 0;
}



この投稿にコメントする

削除パスワード

発言に関する情報 題名 投稿番号 投稿者名 投稿日時
<子記事> Re:プログラムニュートンラフソン法 29234 たかぎ 2006/12/13 17:28:20


No.29234

Re:プログラムニュートンラフソン法
投稿者---たかぎ(2006/12/13 17:28:20)
http://takagi.in/


偶然下記で同じ質問をされている方がいるようです。
参考にされてはいかがでしょうか?
http://okwave.jp/qa2598408.html



この投稿にコメントする

削除パスワード

No.29236

Re:プログラムニュートンラフソン法
投稿者---ヒロ4(2006/12/13 18:33:30)


>偶然下記で同じ質問をされている方がいるようです。
>参考にされてはいかがでしょうか?

完璧私ですね。苦笑


この投稿にコメントする

削除パスワード

No.29318

Re:プログラムニュートンラフソン法
投稿者---かずま(2006/12/25 14:25:45)


http://www.iis.it-hiroshima.ac.jp/~kubokawa/2005/algorithm2/algorithm2-11.pdf
に詳しい説明があるので、その通りにやれば何でもない問題ですね。
もっとも、逆行列を知らなければどうにもなりませんが。
#include <iostream>
#include <iomanip>

double f1(double x1, double x2) { return x1*x1 + x2*x2 - 1; }
double f2(double x1, double x2) { return (x1+x2)*(x1+x2) - 4*x1 - 1; }

double df1x1(double x1, double x2) { return 2 * x1; }
double df1x2(double x1, double x2) { return 2 * x2; }
double df2x1(double x1, double x2) { return 2 * (x1 + x2 - 2); }
double df2x2(double x1, double x2) { return 2 * (x1 + x2); }

void newton_raphson(double x1, double x2)
{
    using namespace std;
    const double eps = 1e-6;
    double dx1, dx2, j;
    cout << fixed;
    for (int i = 0; i < 100; i++, x1 += dx1, x2 += dx2) {
        cout << setw(3) << i << ": x1 =" << setw(10) << x1
             << ", x2 =" << setw(10) << x2 << endl;
        j = df1x1(x1,x2) * df2x2(x1,x2) - df1x2(x1,x2) * df2x1(x1,x2);
        dx1 = (-df2x2(x1,x2) * f1(x1,x2) + df1x2(x1,x2) * f2(x1,x2)) / j;
        dx2 = ( df2x1(x1,x2) * f1(x1,x2) - df1x1(x1,x2) * f2(x1,x2)) / j;
        if (dx1*dx1 + dx2*dx2 < eps * eps) break;
    }
}

int main(void)
{
    newton_raphson(1, 1);
    std::cout << "---\n";
    newton_raphson(1, -1);
    return 0;
}



この投稿にコメントする

削除パスワード

管理者用メニュー    ツリーに戻る    携帯用URL    ホームページ    ログ    タグ一覧