質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

90.42%

  • C

    4110questions

    C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

C言語初心者です。磁場の計算にて、台形公式を用いた重積分がうまくいきません。

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 114

tohokuplug

score 3

永久磁石がある点に作る磁束密度をC言語にて計算したいと考えています。

計算がうまくいかず、どこに間違いがあるのか見つけられていないので、
電磁気ないしC言語に見識のある方のお力添えを頂きたく存じます。

表面磁場 M = 0.24 T 
磁石長さ L = 1 mm (y軸方向)
磁石幅 W = 1 mm (z軸方向)
の磁石がある点に作る磁場を求めます。

磁石表面がy-z表面上にあるとして、長さ方向(y軸方向)にn分割してdy、幅方向(z軸方向)にm分割してdzとし、
それぞれの領域が作る磁場をyとzに関する重積分で重畳して求めようとしています。

この時、磁石の表面磁束密度を M で一様だと仮定すると、分割された各要素の
持つ磁極の強さは dy*dz*M  [Wb] となります。

したがって、磁場のクーロンの法則より、各要素から r [m]離れた点に磁石の要素が作る磁場dHは、

dH = (1/4πµo) * {(dy*dz*M) / r^2}

磁束密度dB [T]に直して、

dB = (1/4π) * {(dy*dz*M) / r^2)}

ここで、
要素の座標が ( xs0 , ys0 , zs0 ) (添え字はソース点の意味でsです。読みにくくてすみません)
磁束密度を求める点が ( x , y , z )
だとすると、2点間の距離 r = √{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2} なので

dB = (1/4π) * [(dy*dz*M) /{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}]

さらに、x軸方向の磁束密度を求めたいため

dB = (1/4π) * [(dy*dz*M) /{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}] * [(x-xs0) / √{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}]

としてそれぞれ求め、これを
y軸方向に 0 → 10^(-3) [m]
z軸方向に 0 → 10^(-3) [m]
の範囲で積分して計算しています。

ここまでで間違いはございますでしょうか。
簡単な静磁場の計算ですが、得手でなくお恥ずかしい限りです。

積分には台形公式を用いています。以下、実際のソースコードになります。
磁束密度を求める点は ( 10^(-3) , 0 , 0 )としています(磁石表面の原点から垂直方向に1 mm の位置)

このプログラムで計算すると、エラーは出ないのですが、解が異常に小さく出たり、
分割数を変えると解の大きさが1桁単位で変化したりと、明らかに間違っており、
おそらくどこか基礎的なことで躓いていることは予想できるのですが、
具体的にどこが間違っているのかが分かりません。

見識ある方のお力添えを頂きたく存じます。

#define _USE_MATH_DEFINES
#include<stdio.h>
#include<math.h>
#define f(ys0,zs0) ((x-xs0)*M)/(4*M_PI)/((x-xs0)*(x-xs0)+(y-ys0)*(y-ys0)+(z-zs0)*(z-zs0))/sqrt((x-xs0)*(x-xs0)+(y-ys0)*(y-ys0)+(z-zs0)*(z-zs0))

int main(void) {
    int n, m, i;
    double L1, L2, W1, W2, M, xs0, xs1, x, ys0, ys1, y, zs0, zs1, z, dy, dz, dSy1, dSy2, Sy1, Sy2, dV, V;

    //要素分割数
    n = 1000;
    m = 1000;
    //-- パラメータの設定
    L1 = 0.000;//yの下限値(磁石の長さ)
    L2 = 0.001;//yの上限値
    W1 = 0.000;//zの下限値(磁石の幅)
    W2 = 0.001;//zの上限値
    M = 0.240;//磁石の磁界強度 [T] = [Wb/m^2]
    xs0 = 0.000;//ソース点
    x = 0.001;//目標点
    y = 0.000;//目標点
    z = 0.000;//目標点
    //--
    //微小量の定義
    dy = (L2 - L1) / n;
    dz = (W2 - W1) / m;
    //初期化
    V = 0;

    //-- 積分
    for (i = 1; i < (n-1); i++) {
        //y軸ソース位置の指定
        ys0 = L1 + dy * i;
        ys1 = ys0 + dy;
        //ループ毎の初期化
        Sy1 = 0;
        Sy2 = 0;
        for (i = 1; i < (m-1); i++) {
            //z軸ソース位置の指定
            zs0 = W1 + dz * i;
            zs1 = zs0 + dz;
            //y軸方向の微小面積計算
            dSy1 = (f(ys0, zs0) + f(ys0, zs1)) * dz / 2;
            dSy2 = (f(ys1, zs0) + f(ys1, zs1)) * dz / 2;
            printf("dSy1 = %lf\ndSy2 = %lf\n", dSy1, dSy2);
            //y軸方向の面積総和
            Sy1 += dSy1;
            Sy2 += dSy2;
        }
        //z軸方向の微小体積積分
        dV = (Sy1 + Sy2)*dy / 2;
        //z軸方向の体積総和
        V += dV;
    }
    //-- 積分終了

    printf("磁束密度 %lf [T]\n", V);
        getchar();
}

試したこと

分割数の変更、各パラメータの値の変更、台形公式の見直し、簡単な関数での重積分結果と比較など

補足情報(FW/ツールのバージョンなど)

Microsoft Visual Studio 2017 を使用

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • KSwordOfHaste

    2019/01/11 23:55 編集

    > ※は見出し表示の回避
    これは正しい解決法ではありません。見出し表示を回避できたとしても字下げが再現されません。
    https://teratail.com/help#about-markdown の「コードを入力」の項に正しい書き方があります。それを参考に質問文を編集ください。質問文の投稿・編集のページ右側にプレビューが表示されていると思いますのでそれで適切な表示になっていることを確認できると思います。
    ---
    また、タイトルの「C言語初心者です」は蛇足なので削除ください。その代わり質問編集画面の左上の初心者マークをクリックし、質問自体に初心者マークを付与してください。

    キャンセル

  • tohokuplug

    2019/01/12 17:31

    ご指摘頂き誠にありがとうございます。ご指摘内容について修正いたしました。以降の投稿では気を付けます。

    キャンセル

回答 1

checkベストアンサー

+2

forループがどっちもi

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2019/01/12 17:32

    解決いたしました。ご指摘頂きまして誠にありがとうございます。非常に初歩的なところでしたので、反省しきりです。

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 90.42%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

同じタグがついた質問を見る

  • C

    4110questions

    C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

  • トップ
  • Cに関する質問
  • C言語初心者です。磁場の計算にて、台形公式を用いた重積分がうまくいきません。