（C言語）ニュートン法のプログラムの為の数値微分

Visual StudioでC言語のプログラムを書く課題が出されたのですが、ニュートン法のプログラムを書くために
関数f(x)=3+ 2x + {2x - cos^2(x) - sin(x)} / {x^2 + cos(x^2) + sin^3(x^3)}
という関数を微分しなければいけません。担当の先生からは、関数が複雑で普通に微分するのは難しいから
f'(x) ≒ {f(x+Δx) - f(x)} / Δx を、Δx=0.0001として数値微分しなさいと指定されています。
どのように計算すればいいのでしょうか？
Δxにそのまま代入するだけではとても計算できそうにないです。

質問者からの補足コメント

• 説明不足なばかりにお手数おかけしてしまって本当に申し訳ありません。
  実は、こんなプログラムを改変して作成しろ、と言われていまして・・・
  文字数がオーバーしてしまったので画像を貼らせていただきます。

  
    補足日時：2015/01/30 16:57

• すみません、画像が元々のサイズより小さく表示されて見づらくなってしまいましたので、プログラムを書かせていただきます。
  #include<stdio.h>
  #include<math.h>
  
  #define MAX_ITERATIONS 50　/* ニュートン法の最大反復回数 */
  
  double function(double x);
  double dFunction(double x);
  double newtonSolver(double initValue, double tolerance);
  
  void main(){
  newtonSolver(4.0, 0.001); /* 1つ目の引数が初期値、2つ目は許容誤差 */
  }

    補足日時：2015/01/30 16:59

• /* ニュートン法による解の計算 */
  /* initValue x の初期値 */
  /* tolerane 許容誤差 */
  double newtonSolver(double initValue, double tolerance){
  double x = initValue; /* 初期値の代入 */
  printf("starting Newton solver...\n");
  printf("step\tx\t\tcorrection\t\n");
  
  for (int counter = 0; counter < MAX_ITERATIONS; ++counter){
  /* 求める解への修正量を求める */
  double correction = -function(x) / dFunction(x);

    補足日時：2015/01/30 17:00

• printf("%3d\t%5f\t%f\n", counter, x, correction);
  
  x += correction;
  
  /* 修正量(correction)の絶対値が許容誤差(tolerance)以下だったら
  解となる x を返して関数を終了する
  ちなみにfabs()は絶対値を求める関数 */
  if (fabs(correction) < tolerance){
  printf("ans = %f\n\n", x);
  return x;
  }
  }
  /* 最大反復回数を超えた場合は、収束していないので
  異常値をあらわすNot a Number (NAN)を返す */
  printf("Newton solver did not converge!\n\n");
  return NAN;
  }

    補足日時：2015/01/30 17:01

• /* 対象の関数 */
  double function(double x){
  return 0;
  }
  
  /* 対象の関数を微分した関数 */
  double dFunction(double x){
  return 0;
  }
  
  
  
  このようなプログラムです。この「対象の関数」と「対象の関数を微分した関数」をなんとかして編集しなければいけないんだろうは思っているのですが、どうすればいいのでしょうか？本当に申し訳ありませんがよろしくお願いします。

    補足日時：2015/01/30 17:07

• おかげさまでプログラムを走らせることができました。本当にありがとうございます。
  ただ、この関数の解は、-0.559629になるらしいのですが、どこが間違っていると思われますか？
  もしかしたら、sin^3(x^3)の部分なのでしょうか？

    補足日時：2015/01/30 20:15

• すみません、プログラミングは大の苦手で、質問の意味が正しく理解できていないかもしれませんが、一応返答させていただきます。
  解がどの場所での微分値なのかとのことですが、恐らく初期値、x=4ではないかと思います。
  あと、powの定義というのはちょっと分からないです。申し訳ありません。

  No.9の回答に寄せられた補足コメントです。 補足日時：2015/01/30 22:53

• 関数f(x)=3+ 2x + {2x - cos^2(x) - sin(x)} / {x^2 + cos(x^2) + sin^3(x^3)}が間違いで、正しいのは
  関数f(x)=3+ 2x + {2x - cos^2(x) + sin(x)} / {x^2 + cos(x^2) + sin^3(x^3)}でした。そのせいで答えが大きく違ってしまったようです。
  あとは/*対象の関数*/
  double function(double x){
  double a = pow(x, 2);
  double b = pow(x, 3);
  return 3 + 2 * x + (2 * x - pow(cos(x), 2) + sin(x)) / (pow(x, 2) + cos(a) + pow(sin(b),3));
  }
  で、解が-0.559629となりました。

    補足日時：2015/01/30 23:29

• お二方にはご迷惑をおかけしました・・・
  これを機に、少しプログラミングの勉強をしようと思います。本当に助かりました。ありがとうございます。

    補足日時：2015/01/30 23:30

No.8 ベストアンサー 回答者： NemurinekoNya 回答日時： 2015/01/30 21:36

#include<stdio.h>

#include<math.h>

//#define MAX_ITERATIONS 50 　/* ニュートン法の最大反復回数 */
#define MAXITERATIONS 50 //上だとエラーになるのでこう書き換えた

double function(double x);
double dFunction(double x);
double newtonSolver(double initValue, double tolerance);

void main(){
double x;
x = newtonSolver(4.0, 0.001); /* 1つ目の引数が初期値、2つ目は許容誤差 */
printf("function(%f) = %f\n", x, function(x)); // ０になるかのチェックだよ
}

/* ニュートン法による解の計算 */
/* initValue x の初期値 */
/* tolerane 許容誤差 */
double newtonSolver(double initValue, double tolerance){
double x = initValue; /* 初期値の代入 */
int counter; // C言語だとこうしないといけない

printf("starting Newton solver...\n");
printf("step\tx\t\tcorrection\t\n");

/* for (int counter = 0; counter < MAXITERATIONS; ++counter){ これはC++やJavaだ */
for (counter = 0; counter < MAXITERATIONS; ++counter){
/* 求める解への修正量を求める */
double correction = -function(x) / dFunction(x);
printf("%3d\t%5f\t%f\n", counter, x, correction);
x += correction;

/* 修正量(correction)の絶対値が許容誤差(tolerance)以下だったら
解となる x を返して関数を終了する
ちなみにfabs()は絶対値を求める関数 */
if (fabs(correction) < tolerance){
printf("ans = %f\n\n", x);
return x;
}
}
/* 最大反復回数を超えた場合は、収束していないので
異常値をあらわすNot a Number (NAN)を返す */
printf("Newton solver did not converge!\n\n");
return NAN;
}

/*対象の関数*/
double function(double x){
return 3 + 2*x + (2*x - pow(cos(x),2) - sin(x)) / (pow(x,2) + cos(pow(x,2)) + pow(pow(sin(x), 3),3));
}

//関数の部分はNo.3さんからのコピー（ポリポリ）
//間違っていたら、No.3さんのせいだからね～（笑）

/* 対象の関数を微分した関数 */
double dFunction(double x){
double dx;
dx = 0.0001;
return (function(x+dx)-function(x))/dx;
}



私が動かしたプログラムはこれ↑よ。
ダメならば、上のソースプログラムをコピペして、コンパイルしろ！！

計算結果は

starting Newton solver...
step x correction
0 4.000000 -6.135696
1 -2.135696 2.833077
2 0.697381 -1.283641
3 -0.586260 -0.494514
4 -1.080774 0.095784
5 -0.984990 0.014808
6 -0.970182 0.000271
ans = -0.969911

function(-0.969911) = -0.000000

です。

No.10 回答者： NemurinekoNya 回答日時： 2015/01/30 23:21

C系のコンピュータ言語では

　a^3 = pow(a,3)
という風に書くんですよ。
　a^n = pow(a,n)
ね。

ただ、文法にうるさいCコンパイラだと、
　pow(a,b)
のa、bのデータの型が実数でなければならないので、bに整数型のデータを入れるとエラーにされてしまう。
BorlandのCコンパイラーなどがその代表。

難しく考えないで、
この質問ページにある

return 3 + 2*x + (2*x - pow(cos(x),2) - sin(x)) / (pow(x,2) + cos(pow(x,2)) + pow(pow(sin(x), 3),3));

をドラッグ＆ドロップし、コピーして、ソースプログラムの問題箇所に貼り付ければいいんですよ。


ドラッグ＆ドロップ
http://www.pc-master.jp/words/drag-drop.html

コピー＆ペースト
http://siriasu.s10.xrea.com/etc/copipe.htm


それで動くと思うんだけれどな。

No.9 回答者： remokon 回答日時： 2015/01/30 22:40

No.1です。

NemurinekoNyaさんが複数回返答されてますね。


＞この関数の解は、-0.559629になるらしいのですが、
は、どの場所での微分値でしょうか。

例えば、x=3での、与えられた条件の微分値は以下です。

f:3+ 2*x + (2*x - cos(x)^2 - sin(x)) / (x^2 + cos(x^2) + sin(x^3)^3);
define(F(x),f);
x1:3$
dx:0.0001$
(F(x1+dx)-F(x1))/dx, numer;

(%i3) f:(-sin(x)-cos(x)^2+2*x)/(sin(x^3)^3+cos(x^2)+x^2)+2*x+3
(%i4) define(F(x),f)
(%i5) x1:3
(%i6) dx:1.0E-4
(%i7) ev((F(dx+x1)-F(x1))/dx,numer)
(%o7) 3.407411997962129

上はMaximaでの実行例です。
MaximaにC形式の数式を出力する機能は、私の知りうる範囲では、ありません。
Fortranならあります。以下です。

(%i4) fortran(f)
(-sin(x)-cos(x)**2+2*x)/(sin(x**3)**3+cos(x**2)+x**2)+2*x+3

上のfが
＞関数f(x)=3+ 2x + {2x - cos^2(x) - sin(x)} / {x^2 + cos(x^2) + sin^3(x^3)}
と一致しているかどうかを確認して下さい。
そして、その-0.559629がどの場所でのxなのか教えて下さい。

なお、pow()はCライブラリmath.h中にあります。貴方のCコンパイラでのpowの定義を教えて下さい。

No.7 回答者： NemurinekoNya 回答日時： 2015/01/30 21:09

☆ただ、この関数の解は、-0.559629になるらしいのですが、どこが間違っていると思われますか？

もしかしたら、sin^3(x^3)の部分なのでしょうか？
◇どこが違いますかといわれましても、
ソースプログラムを見せていただかないことには、
何とも言えませんよ（ニコニコ）。

たぶん、
対象の関数、対象の関数を微分した関数
のどこかを書き間違えているんでしょうね。

ですから、
コピペですよ、コピペ。

return 3 + 2*x + (2*x - pow(cos(x),2) - sin(x)) / (pow(x,2) + cos(pow(x,2)) + pow(pow(sin(x), 3),3));
と
double dx;
dx = 0.0001;
return (function(x+dx)-function(x))/dx;
をコピーして、プログラムに貼り付ける。

何を隠そう、
オレもこの質問文にあるプログラムをコピペしただけだし（ポリポリ）。
質問者さんが補足に書いたプログラムと、おれが書いた部分もコピペした奴をコンパイル、実行しただけだもん（笑）。
そうすると、
ans = -0.969911
になるはずですよ。

この時、
function(-0.969911)がほとんどゼロになることも確かめてあります。


うまくいかなかったら、また連絡してください。

No.6 回答者： NemurinekoNya 回答日時： 2015/01/30 19:09

補足にあったプログラムに、No.5

のプログラムをくっつけコンパイルして、実行させたよ。

私のCコンパイラーはGCCなんで多少の修正をする必要があった。

　───このプログラムのC言語は、CじゃなくC++が混じってる！！───
実行結果は
starting Newton solver...
step x correction
ans = -0.969911

となるようでございます。

No.5 回答者： NemurinekoNya 回答日時： 2015/01/30 18:31

あっ、ゴメン、ゴメン。

NO.4の回答にはウソがあります。


/*対象の関数*/
double function(double x){
return 3 + 2*x + (2*x - pow(cos(x),2) - sin(x)) / (pow(x,2) + cos(pow(x,2)) + pow(pow(sin(x), 3),3));
}

//関数の部分はNo.3さんからのコピー（ポリポリ）
//間違っていたら、No.3さんのせいだからね～（笑）

/* 対象の関数を微分した関数 */
double dFunction(double x){
double dx;
dx = 0.0001;
return (function(x+dx)-function(x))/dx;
}

です。

No.4 回答者： NemurinekoNya 回答日時： 2015/01/30 17:23

/*対象の関数*/

double function(double x){
return 3 + 2*x + (2*x - pow(cos(x),2) - sin(x)) / (pow(x,2) + cos(pow(x,2)) + pow(pow(sin(x), 3),3));
}

//関数の部分はNo.3さんからのコピー（ポリポリ）
//間違っていたら、No.3さんのせいだからね～（笑）


/* 対象の関数を微分した関数 */
double dFunction(double x){
double dx;
dx = 0.0001;
return (f(x+dx)-f(x))/dx;
}


でいいんじゃない。

No.3 回答者： remokon 回答日時： 2015/01/29 20:41

No.2さんの回答を見ると、その場所はプログラム中で指定しても良いようですね。

その場所をx1として、作って見ました。
環境は

Microsoft(R) 32-bit C/C++ Optimizing Compiler Version 15.00.30729.01 for 80x86

です。
これはVisual Studio 2008 Express Editionのコマンドライン版です。


<x.c>

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

float f(double x) {
return x*x;
return 3 + 2*x + (2*x - pow(cos(x),2) - sin(x)) / (pow(x,2) + cos(pow(x,2)) + pow(pow(sin(x), 3),3));
}

void main() {
double x1 = 3;
double dx = 0.00001;

printf("%f\n", (f(x1+dx) - f(x1))/dx);
}


コンパイルは、パスの通った場所で
C:>cl x.c

x.exeを実行すると、次の結果。
6.008148

一般的にはf(x)=x^2として、f'(3)=2x=2*3=6。

No.2 回答者： NemurinekoNya 回答日時： 2015/01/29 18:51

こんばんは。

double f(double x) {
return 問題の関数;
}

とおけば
　df = (f(x+h)-f(x))/h
で数値微分できるでしょっ。
ここでh = Δx = 0.0001

h = 0.0001;
x = 何か適当な値;
while(1) {
df = (f(x+h)-f(x))/h;　　　　//ここが数値微分
xn = x - f(x)/df;
if (収束条件) break;
x = xn; //xの更新
}

で、お求めのニュートン法の本体部分が出来上がり。
dfが0になる場合の例外処理などが必要になりますけれども、これで計算できるはず。
df = f'(x)
ね。

先生は、こういう風に計算をしなさいとおっしゃている。


不安ならば、反復制限を入れる。
whili文の前に
　n = 0;
を入れて、
while (++n < 100)
とか
　while(n++ < 100)
とかすればよい。

このあたりは好みがあるので、while文ではなく
　for(n = 0; n < 100; n++)
でも何でもいいです。

do while文を使ってもいい。

一般的にニュートン法は収束が速いから、10回も反復させれば所定精度の数値解は求まると思いますけれど。

No.1 回答者： remokon 回答日時： 2015/01/29 16:53

そのf(x)が

f:3+ 2*x + (2*x - cos(x)^2 - sin(x)) / (x^2 + cos(x^2) + sin(x^3)^3);

と同じであるなら、計算ソフトで微分できるようです。


(%i3) f:(-sin(x)-cos(x)^2+2*x)/(sin(x^3)^3+cos(x^2)+x^2)+2*x+3
(%i4) diff(f,x)
(%o4) (2*cos(x)*sin(x)-cos(x)+2)/(sin(x^3)^3+cos(x^2)+x^2)
-(-sin(x)-cos(x)^2+2*x)*(9*x^2*cos(x^3)*sin(x^3)^2-2*x*sin(x^2)+2*x)/(sin(x^3)^3+cos(x^2)+x^2)^2+2


よって一般解は上ですが、先生が要求しているのは、特定の場所での微分値を求める作業だと思います。その特定の場所は何処ですか。