掲示板利用宣言

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

 私は

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

掲示板2

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

No.28251

3変数のときのNewton-Raphson法
投稿者---motu(2006/09/23 05:06:51)


Newton-Raphson法を使って3変数の連立方程式を解くプログラムを作ったのですが、うまくいきません。1変数のときはうまくいくのですが…原因が分りません。アドバイスお願いします。

テキストファイルからデータを読み込むようにしてあります。
・3変数のとき

#include <stdio.h> /* インクルードファイル挿入 */
#include <conio.h>
#include <math.h>
#include <stdlib.h>
#define MAXCHR 100 /* 記号定義 c1 */
#define MAXCHR1 101
#define MAXK 100 /* kmax 最大値 */
void main() { /* main 関数開始 */
/* 変数宣言 */
FILE *Finp; /* 入力ファイルポインタ c2 */
FILE *Fout; /* 出力ファイルポインタ */
char cbuf[MAXCHR1]; /* 文字列 */
char Inpfile[]="C:\\Documents and Settings\\哲\\siken2.dat";
/* 入力ファイル名 c3 */
char Outfile[]="C:\\Documents and Settings\\哲\\siken2.out";
/* 出力ファイル名 c4 */
/*----------------------------------------------------------------------*/
int k; /* 整数 汎用インデックス c5 */
int kmax; /* k=1,2,.,kmax(<=MAXK) */
double eps; /* 収束判定イプシロンε c6 */
double a,b,c; /* 2 次方程式パラメータ */
double f1,f2,f3,df11,df12,df13,df21,df22,df23,df31,df32,df33; 
double x0p,y0p,z0p; /* 初期値 */
double x,y,z,xold,yold,zold;
/*----------------------------------------------------------------------*/
printf("プログラム開始、リターンキーを押してください。\n");
getch(); /* 画面へのメッセージ c7 */
printf("\n計算を開始しました。\n");
if ((Finp=fopen(Inpfile,"rt"))==NULL) { /* 入力ファイルオープン c8 */
printf("\n 入力ファイル Example4.dat が開けません。");
exit(1);
}
if ((Fout=fopen(Outfile,"wt"))==NULL) { /* 出力ファイルオープン c9 */
printf("\n 出力ファイル Example4.out が開けません。");
exit(2);
}
/*----------------------------------------------------------------------*/
printf("\n根をNewton-Raphson 法で求めます。");
/* データの読み込み */
fgets(cbuf,MAXCHR,Finp); /* 一行分読み捨て c10 */
fscanf(Finp,"%d%lf%lf%lf%lf%lf%lf%lf",&kmax,&eps,&a,&b,&c,&x0p,&y0p,&z0p);
/* kmax,eps,a,b,c,x0p,y0p,z0p 読込c11 */
if (kmax<1) { kmax=1; } /* 1≦kmax≦MAXK c12 */
if (kmax>MAXK) { kmax=MAXK; }
/*----------------------- a*x^2+b*x+c=0 の解 ----------------------------*/
fprintf(Fout,"Example4 Newton-Raphson 法による解");
fprintf(Fout,"\n\na=%lf b=%lf c=%lf kmax=%d eps=%le xop=%lf yop=%lf zop=%lf",a,b,c,kmax,eps,x0p,y0p,z0p);
fprintf(Fout,"\n\n k x y z"); /* タイトルなど c13 */
k=0;
xold=x0p;
yold=y0p; 
zold=z0p;/* 初期値x0 y0 z0 c14 */
fprintf(Fout,"\n %5d %20.16lf %20.16lf %20.16lf",k,xold,yold,zold);
for (k=1;k<=kmax;k++) { /* k=1,2,...,kmax c15 */
f1=a*xold*xold+b*yold+c*zold*zold+1; /* f(xold) c16 */
df11=2*a*xold; /* f1をxについて偏微分 c17 */
df12=b;
df13=2*c*zold;
f2=b*xold*xold+c*yold+a*zold*zold+2;
df21=2*b*xold;
df22=c;
df23=2*a*zold;
f3=c*xold*xold+a*yold+b*zold*zold+3;

df31=2*c*xold;
df32=a;
df33=2*b*zold;
x=xold-(1/((df11*df22*df33)+(df21*df32*df13)+(df31*df12*df23)-(df11*df32*df23)-(df31*df22*df13)-(df21*df12*df33)))*(((f1*(df22*df33-df23*df32)+f2*(df32*df13-df33*df12)+f3*(df12*df23-df13*df22))));
y=yold-(1/((df11*df22*df33)+(df21*df32*df13)+(df31*df12*df23)-(df11*df32*df23)-(df31*df22*df13)-(df21*df12*df33)))*(((f1*(df23*df31-df21*df33)+f2*(df33*df11-df31*df13)+f3*(df13*df21-df11*df23))));
z=zold-(1/((df11*df22*df33)+(df21*df32*df13)+(df31*df12*df23)-(df11*df32*df23)-(df31*df22*df13)-(df21*df12*df33)))*(((f1*(df21*df32-df22*df31)+f2*(df31*df12-df32*df11)+f3*(df11*df22-df12*df21))));
if (fabs(df11)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df11>=0){
df11=1.0e-20;
} else {
df11=-1.0e-20;
}
}
if (fabs(df12)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df12>=0){
df12=1.0e-20;
} else {
df12=-1.0e-20;
}
}
if (fabs(df13)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df13>=0){
df13=1.0e-20;
} else {
df13=-1.0e-20;
}
}
if (fabs(df21)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df21>=0){
df21=1.0e-20;
} else {
df21=-1.0e-20;
}
}
if (fabs(df22)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df22>=0){
df22=1.0e-20;
} else {
df22=-1.0e-20;
}
}
if (fabs(df23)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df23>=0){
df23=1.0e-20;
} else {
df23=-1.0e-20;
}
}
if (fabs(df31)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df31>=0){
df31=1.0e-20;
} else {
df31=-1.0e-20;
}
}
if (fabs(df32)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df32>=0){
df32=1.0e-20;
} else {
df32=-1.0e-20;
}
}
if (fabs(df33)<1.0e-20) { /* ゼロ割の除去 c18 */
if (df33>=0){
df33=1.0e-20;
} else {
df33=-1.0e-20;
}
}
fprintf(Fout,"\n %5d %20.16lf %20.16lf %20.16lf",k,x,y,z);
if (fabs(x-xold)<eps) { /* 収束条件 c20 */
break;
}
if (fabs(y-yold)<eps) { /* 収束条件 c20 */
break;
}
if (fabs(z-zold)<eps) { /* 収束条件 c20 */
break;
}
xold=x;
yold=y;
zold=z;
}
/*-----------------------------------------------------------------------*/
printf("\n\n 計算が全て終了しました。何かキーを押してください。\n");
if (fclose(Finp)) {
printf("\n 入力ファイルは閉じられませんでした。");
} /* 入力ファイルクローズ */
if (fclose(Fout)) {
printf("\n 出力ファイルは閉じられませんでした。");
} /* 出力ファイルクローズ */
getch(); /* 何かキーが押されるのを待っている。*/
}

データ

Example4 Newton-Raphson 法 2000.5.20
100 1.e-15 1.0 3.0 2.0 -5.0 -5.0 -5.0
kmax eps a b c x0p y0p z0p


長々とすみません。


この投稿にコメントする

削除パスワード

発言に関する情報 題名 投稿番号 投稿者名 投稿日時
<子記事> Re:3変数のときのNewton-Raphson法 28252 επιστημη 2006/09/23 07:26:01
<子記事> Re:3変数のときのNewton-Raphson法 28255 nano 2006/09/23 09:26:58
<子記事> Re:3変数のときのNewton-Raphson法 28258 nano 2006/09/23 13:26:12


No.28252

Re:3変数のときのNewton-Raphson法
投稿者---επιστημη(2006/09/23 07:26:01)


>Newton-Raphson法を使って3変数の連立方程式を解くプログラムを作ったのですが、うまくいきません。1変数のときはうまくいくのですが…原因が分りません。アドバイスお願いします。

インデントされてない、なにがどううまくいかんか説明なし、長すぎ。
「読んでくれなくて結構です」と言っているに等しい。



この投稿にコメントする

削除パスワード

No.28255

Re:3変数のときのNewton-Raphson法
投稿者---nano(2006/09/23 09:26:58)


>Newton-Raphson法を使って3変数の連立方程式を解くプログラム

ソースの中のコメントでは2次方程式を解くようなことが
書いてありますが、結局何がしたいのでしょうか?



この投稿にコメントする

削除パスワード

No.28258

Re:3変数のときのNewton-Raphson法
投稿者---nano(2006/09/23 13:26:12)


もし、x^2 + 3*x + 2 = 0 という2次方程式の解の1つを求めたかったのであれば、
例えばこんな感じでしょうか。
#include <stdio.h>
#include <math.h>

#define EPS       (1.0e-05)
#define ITERATION (100)

double f(double x);
double g(double x);

int main(void)
{
    double x0 = 0.0, x1, g0;
    int i;
    
    for (i = 0; i < ITERATION; i++) {
        g0 = (g(x0) == 0.0) ? EPS : g(x0);
        x1 = x0 - f(x0) / g0;
        if (fabs(x1 - x0) < EPS)
            break;
        x0 = x1;
    }
    if (i == ITERATION)
        printf("収束しませんでした。\n");
    else 
        printf("解=%.6f\n", x1);
    return 0;
}

double f(double x)
{
    return x * x + 3 * x + 2;
}

double g(double x)
{
    return 2 * x + 3;
}




この投稿にコメントする

削除パスワード

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