ここから本文です

! n 個の質点からなる連成振動子の運動を 0 < t < tend の範囲で解く。長さ 2n ...

アバター

ID非公開さん

2019/7/2014:20:10

! n 個の質点からなる連成振動子の運動を 0 < t < tend の範囲で解く。長さ 2n の配列 w(1:n, 1:2) を定義して, ! j 番目の質点の変位を w(j, 1),速度を w(j, 2) に格納する。質点の質量の逆数に比例する無次元

パラメータは
! 配列 c(1:n) に格納する。初期条件は変位ゼロ,速度は 1 番目の質点のみ vzero とする。
program coupled_oscillators
implicit none
integer(4), parameter :: n = 10
integer(4) :: j
real(8) :: t, dt
real(8), dimension(n) :: c
real(8), dimension(n, 2) :: w, p, k1, k2, k3, k4 real(8), parameter :: vzero = 0.1d0, tend = 100.0d0 ! 時間ステップの刻み幅 dt を標準入力から読み込む。 read(*, *) dt
! 質量はみんな同じ。
c(:) = 1.0d0
! 未知変数 w と時刻 t の初期化。
w(:, 1) = 0.0d0
w(1, 2) = vzero
w(2:n, 2) = 0.0d0
t = 0.0d0
! 時間積分の開始(無限 do loop)。
do
! 時刻 t での質点の位置を書き出す。
write(*, "(f12.7, 99e16.8)") t, (w(j, 1) + dble(j)/dble(n+1), j = 1, n)
! 時刻が tend を超えたら無限 do loop から脱出する。
if ( t >= tend ) stop
! 4次ルンゲ・クッタ法。微分方程式を dw/dt = f(w) だと思って,機械的にアルゴリズムを適用する。 ! まず補助変数 k1 = f(p) を計算する。
p(:, :) = w(:, :)
call rhs (k1, p, c, n)
! 補助変数 k2 = f(p) を計算する。


これをc言語に変換してください。お願いします。

閲覧数:
30
回答数:
1
お礼:
25枚

違反報告

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

プロフィール画像

カテゴリマスター

uso8megaさん

編集あり2019/7/2015:07:07

/* 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
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/(n+1.0));
    /* 時刻が 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) を計算する。*/

/****残念!尻切れトンボ  後は自力でどうぞ****/

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

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

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

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

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

閉じる

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

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

閉じる