ここから本文です

! 補助変数 k2 = f(p) を計算する。

アバター

ID非公開さん

2019/7/2418:32:34

! 補助変数 k2 = f(p) を計算する。

p(:, :) = w(:, :) + 0.5d0 * dt * k1(:, :)
call rhs (k2, p, c, n)
! 補助変数 k3 = f(p) を計算する。
p(:, :) = w(:, :) + 0.5d0 * dt * k2(:, :)
call rhs (k3, p, c, n)
! 補助変数 k4 = f(p) を計算する。
p(:, :) = w(:, :) + dt * k3(:, :)
call rhs (k4, p, c, n)
! 未知変数 w と時刻 t を次の時間ステップでの値で上書きする。
w(:, :) = w(:, :) + (dt/6.0d0) * (k1(:, :) + k4(:, :) + 2.d0*(k2(:, :) + k3(:, :)))
t = t + dt
end do end
! 微分方程式 dw/dt = f(w) の右辺を計算するサブルーチン。 subroutine rhs (f, w, c, n)
implicit none
integer(4) :: n, j
real(8), dimension(n) :: c
real(8), dimension(n, 2) :: f, w
do j = 1, n
! duj/dt = vj
f(j, 1) = w(j, 2)
! dvj/dt = −cj (2uj − uj−1 − uj+1) if ( j == 1 ) then
f(j, 2) = - c(j) * (2.d0*w(j, 1) - w(j+1, 1))
else if ( j == n ) then
f(j, 2) = - c(j) * (2.d0*w(j, 1) - w(j-1, 1))
else
f(j, 2) = - c(j) * (2.d0*w(j, 1) - w(j-1, 1) - w(j+1, 1))
end if
end do
return end

これをc言語で表して下さい。お願いします。

閲覧数:
22
回答数:
3

違反報告

ベストアンサーに選ばれた回答

プロフィール画像

カテゴリマスター

uso8megaさん

2019/7/2621:32:06

https://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q102107975...
↑の続きですね。ちなみに
↑は少々誤りがありましたの
で今回のに差し替えて下さい。
http://ideone.com/ZXQY49#Fortran元ネタ

/*http://ideone.com/1S7r4O#C変換版*/
/* PROGRAM coupled_oscillators */
/* n 個の質点からなる連成振動子の運動を         */
/* 0 < t < TEND の範囲で解く。                  */
/* 長さ 2*N の配列 w[0~(N-1)][0~1] を定義して, */
/* j 番目の質点の変位を w[j][0],               */
/* 速度を w[j][1] に格納する。                  */
/* 質点の質量の逆数に比例する無次元パラメータは */
/* 配列 c[0~(N-1)] に格納する。                 */
/* 初期条件は変位ゼロ,                         */
/* 速度は 1 番目の質点のみ VZERO とする。       */
#include<stdio.h>
#define N      10
#define VZERO   0.1
#define TEND  100.0
void rhs(double f[N][2], double w[N][2], double c[N], int n);
int main(void){
  int j;
  double t,dt, c[N], w[N][2],p[N][2],
    k1[N][2], k2[N][2], k3[N][2], k4[N][2];
  /* 時間ステップの刻み幅 dt を標準入力から読み込む。*/
  scanf("%lf",&dt);
  /* 質量はみんな同じ。*/
  for(j=0;j<N;++j)c[j] = 1.0;
  /* 未知変数 w と時刻 t の初期化。*/
  for(j=0;j<N;++j) w[j][0] = 0.0,
    w[j][1] = j==0 ? VZERO : 0.0;
  t = 0.0;
  /* 時間積分の開始(無限 for loop)。*/
  for(;;){
    /* 時刻 t での質点の位置を書き出す。*/
    printf("%12.7f",t);
    for(j=0;j<N;++j)
      printf("%16.8e", w[j][0] + (j+1.0)/(N+1.0));
    printf("\n");
    /* 時刻が TEND を超えたら無限 for loop から脱出する。*/
    if ( t >= TEND ) break;
    /* 4次ルンゲ・クッタ法。微分方程式を */
    /* dw/dt = f(w) だと思って,機械的に */
    /* アルゴリズムを適用する。*/
    /* まず補助変数 k1 = f(p) を計算する。*/
    for(j=0;j<N;++j)
      p[j][0]=w[j][0],
      p[j][1]=w[j][1];
    rhs(k1, p, c, N);
    /* 補助変数 k2 = f(p) を計算する。*/
    for(j=0;j<N;++j)
      p[j][0]=w[j][0] + 0.5*dt*k1[j][0],
      p[j][1]=w[j][1] + 0.5*dt*k1[j][1];
    rhs(k2, p, c, N);
    /* 補助変数 k3 = f(p) を計算する。*/
    for(j=0;j<N;++j)
      p[j][0]=w[j][0] + 0.5*dt*k2[j][0],
      p[j][1]=w[j][1] + 0.5*dt*k2[j][1];
    rhs(k3, p, c, N);
    /* 補助変数 k4 = f(p) を計算する。*/
    for(j=0;j<N;++j)
      p[j][0]=w[j][0] + dt*k3[j][0],
      p[j][1]=w[j][1] + dt*k3[j][1];
    rhs(k4, p, c, N);
    /* 未知変数 w と時刻 t を次の時間ステップでの値で上書きする。*/
    for(j=0;j<N;++j)
      w[j][0]=w[j][0] + (dt/6.0)*(k1[j][0]+k4[j][0]+2.0*(k2[j][0]+k3[j][0])),
      w[j][1]=w[j][1] + (dt/6.0)*(k1[j][1]+k4[j][1]+2.0*(k2[j][1]+k3[j][1]));
    t = t + dt;
  }
return 0;}

/* 微分方程式 dw/dt = f(w) の右辺を計算するサブルーチン。*/
void rhs(double f[N][2], double w[N][2], double c[N], int n){
  int j;
  for(j=0;j<n;++j){
    /* du(j)/dt = v(j) */
    f[j][0] = w[j][1];
    /* dv(j)/dt = -c(j) * (2u(j) - u(j-1) - u(j+1)) */
    if ( j == 0 ) {
      f[j][1] = - c[j] * (2.0*w[j][0] - w[j+1][0]);
    } else if ( j == n-1 ) {
      f[j][1] = - c[j] * (2.0*w[j][0] - w[j-1][0]);
    } else {
      f[j][1] = - c[j] * (2.0*w[j][0] - w[j-1][0] - w[j+1][0]);
    }
  }
}

アバター

質問した人からのコメント

2019/7/27 02:11:17

とても丁寧にありがとうございます。
助かりました!

ベストアンサー以外の回答

1〜2件/2件中

並び替え:回答日時の
新しい順
|古い順

2019/7/2618:31:03

Fortran ぽいけどでも私わかりませ~ん

meigr01さん

2019/7/2418:47:10

表示しているコードの言語を示してください。

返信を取り消しますが
よろしいですか?

  • 取り消す
  • キャンセル

この質問につけられたタグ

みんなで作る知恵袋 悩みや疑問、なんでも気軽にきいちゃおう!

Q&Aをキーワードで検索:

Yahoo! JAPANは、回答に記載された内容の信ぴょう性、正確性を保証しておりません。
お客様自身の責任と判断で、ご利用ください。
本文はここまでです このページの先頭へ

「追加する」ボタンを押してください。

閉じる

※知恵コレクションに追加された質問は選択されたID/ニックネームのMy知恵袋で確認できます。

不適切な投稿でないことを報告しました。

閉じる