• 締切済み

C++言語のプログラムをfortranに変換

const int N = 100; const double q = 10.0, dt = 0.00001, Dm = 30.0, t0 = 2.0, K = 1.0, pi = 3.141592, f = 3.0; double C[N], dC[N]; double dx = q/N; for (int i = 0; i < N; ++i) C[i] = 0;　// 初期条件 for (double t = 0; t < t0; t += dt) { 　C[0] = 0;　// 境界条件1 　C[N - 1] = K*sin(2*pi*f*t);　// 境界条件2 　for (int i = 1; i < N - 1; ++i) dC[i] = (Dm*(C[i + 1] - 2*C[i] + C[i - 1])/(dx*dx))*dt; 　for (int i = 1; i < N - 1; ++i) C[i] += dC[i]; } for (int i = 0; i < N; ++i) cout << C[i] << endl; // t = t0 このプログラムをfortranに変換できる方いますか？

みんなの回答

  • FT56F001
  • ベストアンサー率59% (355/599)
回答No.2

だいたい，このくらいです。 処理系によっては，配列の添字は1以上で，0はダメかもしれません。 INTEGER N,IT DOUBLE Q,dt,DM,T0,K,PI,F,dx DATA N/100/,Q/10/,dt/0.00001/,DM/30/,T0/2.0/,K/1.0/,pi/3.141592,/F/3.0/ DOUBLE C(100),DC(100) dx=Q/N * 初期条件 DO 10 I=1,N-1 C(I)=0 10 CONTINUE DO 20 IT=1,T0/dt-1 t=IT*dt C(0)=0 * 境界条件 C(N-1)=K*sin(2.0*pi*f*t) * 境界条件2 DO 30 I=1,N-2 dC(I) = (Dm*(C(I + 1) - 2*C(I) + C(I - 1))/(dx*dx))*dt 30 CONTINUE DO 40 I=1,N-2 C(I)=C(I)+ dC(I) 40 CONTINUE 20 CONTINUE DO 60 I=0,N-1 write(*,*) C(I) 60 CONTINUE * t=T0

te107484p
質問者

お礼

返信ありがとうございます！！！

te107484p
質問者

補足

もう。。。なんて言ったらよいか。。。本当にありがとうございます。 これを基本ベースにして様々な応用をしていきたいと思います！！！

  • joqr
  • ベストアンサー率18% (742/4026)
回答No.1

はい！(元気よく手を挙げる)

te107484p
質問者

補足

申し訳ないのですが。。。お時間のある時にでもしていただきたいのですが。。。 プログラミング詳しいですね(>_<)

  1. カテゴリ
  2. 学問・教育
  3. 自然科学
  4. 物理学

関連するQ&A

  • c言語でDFTのプログラムを作成したのですが

    c言語でDFTのプログラムを作成しました。 以下にソースを載せます。 #include<stdio.h> #include<math.h> #include<stdlib.h> #include<time.h> #define PI 3.141592653589793 #define N 64 //データ数 DFT(double result[]){ int i,k; double A[N],B[N],T=0; //A[N]:実数部,B[N]:虚数部 double a,b; for(k=0;k<N;++k){ a=b=0; for(i=0;i<N;++i){ a=a+result[i]*cos(2*PI*i*k/N);　　　　　 b=b+(-1.0)*result[i]*sin(2*PI*i*k/N); } A[k]=a/N; B[k]=b/N; } for(i=0;i<N;++i){ printf("[%f秒]:Re:%f,Im:%f\n",T,A[i],B[i]); //変換後の値を表示 T=T+(0.1/N); } } main(){ int i; double T=0; double Sampdata; double result[N]; for(i=0;i<N;++i){ Sampdata=5*sin(20*PI*T);　　　　　　//0~0.1秒間をN個にサンプリング result[i]=Sampdata; //サンプリングデータを代入 T=T+(0.1/N); } clock_t start,end; //処理時間計測開始 start=clock(); DFT(result); end=clock(); printf("%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC); //処理時間表示 } 元信号には5sin(20πt)の値を入れています。この信号は周期は0.1secです。 これでフーリエ変換を行うとデータ数N/2を中心に対称なデータが出てくるのですが、処理が終わるのが早い気がするんです。 例えば2^15個のデータで実行しても2分もかからずに処理が終わってしまいます。一応、対称性が出てるとはいえ、終わるのが早すぎる気がするのですが、おかしい所があれば教えていただけると嬉しいです。 よろしくお願いします。

  • フーリエ変換のC言語プログラムについて

    正弦波（およびガウス性雑音)をフーリエ変換(離散)→逆フーリエ変換するというプログラムを組みました。正弦波をフーリエ変換すると実部は２回ピークがくるはずなのですが、すべて「0.000000」または「-0.000000」と表示されてしまいます。虚部は正常なようで実装の仕方もさほど違わないので、何が問題なのかわからずにいます。念のためコードはすべて載せますが、該当箇所は関数Fourierの fp = fopen("reX.txt", "w"); //書き込む あたりです。問題点を教えていただけないでしょうか。お願いします。 //gauss.txt, sin.txt:発生させたガウス性雑音＆正弦波 //reX, imX:フーリエ変換の実部と虚部 //re-1, im-1:逆フーリエ変換の実部と虚部 #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> #define PI 3.14159265358979323846 #define N 256 //n:長さ, w:角周波数, p:位相(phase), a:振幅 void SinCurve(int n, double w, double p, double a) { FILE *fp; double x; int t; fp = fopen("sin.txt", "w"); //書き込むので"w" if(fp == NULL) { printf("file open error\n"); } else { for(t = 0; t < n; t++) { x = a * sin( w*(double)t + p ); fprintf(fp, "%f\n", x); } } fclose(fp); } //n:長さ, s:分散, m:平均 void Gauss(int n, double s, double m) { FILE *fp; double x, x1, x2, y1; int t; srand((unsigned) time(NULL)); fp = fopen("gauss.txt", "w"); //書き込むので"w" if(fp == NULL) { printf("file open error\n"); } else { for(t = 0; t < n; t++) { x1 = ( (double)rand() + 1.0 ) / ( (double)RAND_MAX + 2.0); x2 = ( (double)rand() + 1.0 ) / ( (double)RAND_MAX + 2.0); y1 = pow(-2.0*log(x1), 0.5) * cos(2.0*PI*x2); y1 = s * y1 + m; fprintf(fp, "%f\n", y1); } } fclose(fp); } //ファイル名sのデータをフーリエ変換し、実部と虚部をreX, imXに保存 void Fourier(int num, char *s) { FILE *fp; int k, n; double largeX, x[N+100], t; fp = fopen(s, "r"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { // printf("%s\n", s); for(k = 0; k < num; k++) { fscanf(fp, "%lf", &x[k]); printf("x[%d]=%f\n", k, x[k]); } } fp = fopen("reX.txt", "w"); //書き込む if(fp == NULL) { printf("file open error\n"); } else { for(k = 0; k < num; k++) { largeX = 0.0; t = 2.0*PI*(double)k / (double)N; for(n = 0; n < num; n++) { largeX += x[n] * cos((double)n*t); // printf("%f\n", largeX); } fprintf(fp, "%f\n", largeX); printf("reX[%d]=%f\n", k, largeX); } } fp = fopen("imX.txt", "w"); //書き込む if(fp == NULL) { printf("file open error\n"); } else { for(k = 0; k < num; k++) { largeX = 0.0; t = 2.0*PI*k / (double)N; for(n = 0; n < num; n++) { largeX -= x[n] * sin(n*t); } fprintf(fp, "%f\n", largeX); } } fclose(fp); } void InverseFourier(int num) { FILE *fp; int k, n; double a[N+100], b[N+100], x, t; //a:reX, b:imX fp = fopen("reX.txt", "r"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { for(k = 0; k < num; k++) { fscanf(fp, "%lf", &a[k]); // printf("a[%d]=%f\n", k, a[k]); } } fp = fopen("imX.txt", "r"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { for(k = 0; k < num; k++) { fscanf(fp, "%lf", &b[k]); // printf("b[%d]=%f\n", k, b[k]); } } fp = fopen("re-1.txt", "w"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { for(n = 0; n < num; n++) { x = 0.0; t = 2.0*PI*(double)n / (double)N; for(k = 0; k < num; k++) { x +=a[k] *cos(k*t) - b[k] *sin(k*t); } x /= (double)N; fprintf(fp, "%f\n", x); // printf("x[%d]=%f\n", n, x); } } /* fp = fopen("im-1.txt", "w"); //読み込み if(fp == NULL) { printf("file open error\n"); } else { for(n = 0; n < num; n++) { x = 0.0; for(k = 0; k < num; k++) { t = 2.0*PI*(double)k / (double)N; x = x + a[k] *sin(n*t) + b[k] *cos(n*t); } x /= (double)N; fprintf(fp, "%f\n", x); } } */ fclose(fp); } int main(void) { SinCurve(N, PI/8.0, 0.0, 1.0); // Gauss(N, 1.0, 0.0); Fourier(N, "sin.txt"); // Fourier(N, "gauss.txt"); InverseFourier(N); return 0; }

  • Visual C++での数値計算のプログラミング

    質問初めてになります。 プログラミングにあまり詳しくない大学院の数理科の学生です。 学校で熱方程式の陽解法のプログラミングのレポートが出されたのですが、困っています。 レポートで詰まっている点は windows Visual C++ Expressionでの陽解法のプログラミングです。 indows Visual C++ Expression2008が良くわからないのでコマンドプロンプトで実行しています。 熱方程式の初期値問題の陽解法のプログラムをなんとか組んでいます。 私が組んだプログラムではLinuxでは通るのですがwindows Visual C++ Expressionでは通りません。 このプログラムをwindows Visual C++ Expression　2008で通すにはどのように直したがいいのか教えて頂きたいと思います。よろしくお願いします。 以下は私なりに組んだプログラミングです。 πの値、u1[N+1]、u2[N+2] の３箇所でエラーがでてると思われます。 #include<stdio.h> #include<stdio.h> #include<math.h> int main(){ int i,k,kmax; int N=10; double dx double dt=0.01; double u1[N+1]; double u2[N+1]; double r=dt/(dx*dx); double T=1.0; kmax=T/dt; u1[0]=0; u1[N]=0; for(i=1;i<=N;i++){ 　　u1[i]=sin(M_PI*i*dx); } for(k=1;k<=kmax;k++){ for(i=1;i<=N-1;i++){ u2[i]=r*u1[i-1]+(1-2*r)*u1[i]+r*u1[i+1]; } 　　for(i=1;i<=N-1;i++){ u1[i]=u2[i]; 　　} 　　for(i=0;i<=N;i++){ printf("%f %f %f\n",dt*k,dx*i,u1[i]); 　　} 　　printf("\n"); } return 0; }

  • Ｃ言語のプログラムについてです。

    円を描くプログラムを作りたいのですが、条件としてpointを使わないでつくらなければなりません。あとそのプログラムについての説明もいるので、合わせて回答お願いします。pointを使うと↓のプログラムです。 #include<stdio.h> #include<math.h> int main(void); int main(void) { double point[2][2]; double r; double pi=3.14159; int i; onenpl(); space(-100.0,-100.0,100.0,100.0); printf("円を描きます\n"); printf("半径を入力してください\n"); scanf("%lf",&r); point[0][0]=r*cos(0.0); point[0][1]=r*sin(0.0); for (i=1;i<360;i=i+1) { point[1][0]=r*cos(pi/180.0*i); point[1][1]=r*sin(pi/180.0*i); line(point[0][0],point[0][1],point[1][0],point[1][1]); point[0][0]=point[1][0]; point[0][1]=point[1][1]; } closepl(); return(0); }

  • C++で分からないプログラムがあるんですが

    #include <iostream> #include <cmath> using namespace std; int main() { static const int N = 2; double va[N]={3,-4}; double vb[N]={4,3}; double a,b; double p; for (int i = 0; i < N; ++i) { for (int i = 0; i < N; ++i) { } } cout << "va + vb = (" ; for (int i = 0; i < N; ++i) { cout << va[i] + vb[i]; if (i < N - 1) { cout << ", "; } } cout << ")" << '\n'; cout << "va - vb = (" ; for (int i = 0; i < N; ++i) { cout << va[i] - vb[i]; if (i < N - 1) { cout << ", "; } } cout << ")" << '\n'; p = 0; for (int i = 0; i < N; ++i) { p += va[i] * vb[i]; } cout << "va・vb = " << p << '\n'; a = 0; for (int i = 0; i < N; ++i) { a += va[i] * va[i]; } a = sqrt(a); b = 0; for (int i = 0; i < N; ++i) { b += vb[i] * vb[i]; } b = sqrt(b); if (a * b != 0) { cout << "cosθ = " << p / (a * b) << '\n'; } return 0; } これで、ベクトルの加減とベクトルの内積とcosθが出るんですが、2つのベクトルを適当に初期化しないといけないんですが、初期化ってこれで初期化ってできてますか？

  • C言語のプログラムについて。

    C言語のプログラミングについて質問です。 入力されたデータの配列とデータ数を渡すと配列に格納された値を逆順にして、格納し直す関数reverse関数を書き結果を出力せよ、というものなのですが下のように書いたのですが、うまく作動しません。どこがいけないのでしょうか...？教えていただきたいです。 #include <stdio.h> void reverse(int *data[], int n); #define MAX 100 int main() { int data[MAX]; int n, i; scanf("%d", &n); if (n >= MAX) n = MAX; for (i = 0; i < n; i ++){ scanf("%d", &data[i]); } reverse(data, n); for (i = 0; i < n; i ++) { printf("%d\n", data[i]); } return 0; } void reverse(int *data[], int n) { int c, i; for (i = 0; i < n; i ++) { c = *data[i]; *data[i] = *data[n - (i + 1)]; *data[n - (i + 1)] = c; } }

  • Ｃ言語のプログラムで質問です。

    Ｃ言語のプログラムで質問です。 ヤコビ法で固有値・固有ベクトル絵を求めたいのですが、固有値は上手く出るのですが、 固有ベクトルがめちゃくちゃな値が出ます。何が間違っているのでしょうか教えてください。下のは一応自分で作ってみたプログラムです。 #include <stdio.h> #include <math.h> #define N 4 #define PI 3.14159 void jacobi(double a[N][N], double λ[N], double v[N][N]); int main(void) { int i, j; double a[N][N], λ[N], v[N][N], x; for(i=0;i<N;i++){ printf("\n a[%d][0] a[%d][1] a[%d][2] a[%d][3]=",i,i,i,i); scanf("%lf %lf %lf %lf",&a[i][0],&a[i][1],&a[i][2],&a[i][3]); } jacobi(a, λ, v); for(j=0; j<N; j++){ printf("\n固有値λ[%d]=%lf", j, λ[j]); printf("\n固有ベクトル\n"); for(i=0; i<N; i++){ printf("%lf\n", v[i][j]); } } return 0; } /* ヤコビ法による固有値計算 */ void jacobi(double a[N][N], double λ[N], double v[N][N]) { int i, j, kmax=100, repeat, p, q; double eps, c, s,theta, gmax; double apq, app, aqq, apqmax, apj, aqj, vip, viq; gmax=0.0; for(i=0; i<N; i++){ s=0.0; for(j=i+1; j<N; j++){ s += fabs(a[i][j]); } if(s>gmax) gmax=s; } eps=0.000001*gmax; for(i=0; i<N; i++){ for(j=0; j<N; j++){ v[i][j]=0.0; } v[i][i]=1.0; } for(repeat=1; repeat<kmax; repeat++){ /* 収束判定 */ apqmax = 0.0; for(p=0; p<N; p++){ for(q=0; q<N; q++){ if(p!=q){ apq=fabs(a[p][q]); if(apq>apqmax) apqmax=apq; } } } if(apqmax<eps) break; for(p=0; p<N-1; p++){ for(q=p+1; q<N; q++){ apq=a[p][q]; app=a[p][p]; aqq=a[q][q]; if(fabs(apqmax)<eps) break; /*回転角計算*/ if(fabs(app-aqq)>=1.0e-15){ theta= 0.5*atan(2.0*apq/(app-aqq)); }else{ theta= PI/4.0; } c = cos(theta); s = sin(theta); a[p][p] = app*c*c + 2.0*apq*c*s + aqq*s*s; a[q][q] = app*s*s - 2.0*apq*c*s + aqq*c*c; a[p][q] = 0.0; a[q][p] = 0.0; for(j=0; j<N; j++){ if(j!=p && j!=q){ apj = a[p][j]; aqj = a[q][j]; a[p][j] = apj*c + aqj*s; a[q][j] = -apj*s + aqj*c; a[j][p] = a[p][j]; a[j][q] = a[q][j]; } } /*固有ベクトル*/ for(i=0; i<N; i++){ v[i][p] = v[i][p]*c+v[i][q]*s; v[i][q] = -v[i][p]*s+v[i][q]*c; } } } eps=eps*1.05; } for(i=0; i<N; i++) λ[i] = a[i][i]; }

  • Cのプログラムに変換

    どなたかこのプログラムをCに変化してくれる方いませんか？ なんで変換するのとかの質問はいらないです。 ただ純粋に変換してくれる方いませんか？ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <complex.h> #define NREPEAT 1000 int sweepout(float complex **a, int size) { int n,m,k; // 前進消去 for(n=0; n<size; n++) { // 軸の選択 { float amax = cabs(a[n][n]);; float aabs; int npivot = n; for(m=n+1; m<size; m++) { if((aabs = cabs(a[m][n]))>amax) { npivot = m; amax = aabs; } } if(aabs == 0.0) return n; if(npivot != n) { // 軸の交換 float complex t; for(k=n; k<=size; k++) { t = a[n][k]; a[n][k] = a[npivot][k]; a[npivot][k] = t; } } } // 消去 { float complex t; t = 1.0/a[n][n]; // 対角要素 a[n][n] <-- 1 for(k=n+1; k<=size; k++) a[n][k] *= t; // 非対角要素 a[m][n] <-- 0 for(m=n+1; m<size; m++) { for(k=n+1; k<=size; k++) a[m][k] -= a[m][n]*a[n][k]; } } } // 後退代入 { float t; for(n=size-1; n>0; n--) { for(m=0; m<n; m++) a[m][size] -= a[n][size]*a[m][n]; } } return -1; } // 乱数を発生する double drand() { double rn; rn = (2*((double)random()))/((double)(RAND_MAX))-1.0; return rn; } main(int argc, char **argv) { int n, m, k, ncount; double complex **a0; float complex **a; double complex t; double d, r; if(argc>1) n = atoi(argv[1]); if(n <=0) n = 3; srandom((argc>2)?atoi(argv[2]):0); a0 = (double complex **)malloc(n*sizeof(double complex *)); for(m=0; m<n; m++) a0[m] = (double complex *)malloc((n+1)*sizeof(double complex)); a = (float complex **)malloc(n*sizeof(float complex *)); for(m=0; m<n; m++) a[m] = (float complex *)malloc((n+1)*sizeof(float complex)); for(ncount=0; ncount<NREPEAT; ncount++) { // 乱数により、配列の係数を決める(double型) for(m=0; m<n; m++) for(k=0; k<=n; k++) a0[m][k] = drand() + 1.0I*drand(); // 係数をfloat型に変換し、配列 a に代入する for(m=0; m<n; m++) { for(k=0; k<n; k++) a[m][k] = (float complex)a0[m][k]; t = 0; for(k=0; k<n; k++) t += a0[m][k]*a0[k][n]; a[m][n] = (float complex)t; } // 連立方程式を解き、数値解を a[][n] に求める sweepout(a,n); // 2乗誤差を求める d = r = 0; for(m=0; m<n; m++) { d += pow(cabs(a[m][n]-a0[m][n]),2); r += pow(cabs(a0[m][n]),2); } printf("%lf\t\n", 0.5*log10(d/r)); } for(m=0; m<n; m++) { free(a0[m]); free(a[m]); } free(a0); free(a); }

  • c言語です。

    c言語です。 実行結果 式 3 X1 + 2 X2 + 1 X3 = &g 2 X1 + 5 X2 + 2 X3 = &g 1 X1 + 4 X2 + 1 X3 = &g 解 X1 = 1 X2 = 2 X3 = 3 を 式 3 X1 + 2 X2 + 1 X3 = 10 2 X1 + 5 X2 + 2 X3 = 18 1 X1 + 4 X2 + 1 X3 = 12 解 X1 = 1 X2 = 2 X3 = 3 に直したいのですが&gの所をどのようにしたら10.18.12になりますか？ #include <stdio.h> #include <float.h> #define N 3 double A[N][N] = {{3,2,1}, {2,5,2}, {1,4,1}}; double b[N] = { 10, 18, 12 }; void Gauss_J( int, double*, double* ); void main(void) { int i; printf( "%d式\n", N ); for( i = 0; i < N ; i++ ) { printf( "%g X1 + %g X2 + %g X3 = &g \n", A[i][0], A[i][1], A[i][2], b[i] ); } printf("解\n"); Gauss_J(N, (double *)A, (double *)b ); printf("X1 = %g \n", b[0]); printf("X2 = %g \n", b[1]); printf("X3 = %g \n", b[2]); } void Gauss_J(int n, double *a, double *b) { int p, i, j,I ; double pivot, c ; for ( p = 0 ; p < n ; p++ ) { pivot = a[ p*n + p ]; for ( i = p ; i < n ; i++ ) { a[ p*n + i ] /= pivot; } b[ p ] /= pivot; for ( I = 0 ; I < n ; I++) { if (I != p) { c = a[ I*n + p]; for ( j = p ; j < n; j++ ) { a[ I*n + j] -= c * a[ p*n + j ]; } b[ I ] -= c * b[ p ]; } } } return ; }

  • C言語　正弦関数の級数展開

    はじめまして。 sinxの級数展開を７項まで取った場合と組み込み関数で求めたsinxの値の差を，0度から360度まで，プログラムを作成して求めよ、という問題があるのですが、300度まででとまってしまいます。 どなたか、わかる方がいまいたら、教えてください。 以下リストです。 #include <stdio.h> #include<conio.h> #include <math.h> const double pi=atan(1.0)*4; double kaijyo(int n) { double ret=1; for (int i=2;i<=n;i++) { ret*=i; } return ret; } double polSin(double theta,int order) { int i; double ret=theta; for (i=1;i<order;i++) { if (i%2==1) { ret-=pow(theta,(i*2+1))/kaijyo(i*2+1); } else { ret+=pow(theta,(i*2+1))/kaijyo(i*2+1); } } return ret; } int main() { int order=7; double rad,theta,trueVal,polVal; for (theta=0;theta<360;theta+=1.0) { rad=theta/180*pi; trueVal=sin(rad); polVal=polSin(rad,order); printf("%e\n",polVal-trueVal); } getch(); return 0; }

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

あなたにピッタリな商品が見つかる！ OKWAVEセレクト

質問する