• 締切済み

大学のC言語の課題で

ルンゲクッタ法を用いたバネの単振動を解析する課題をやっています。以下の部分まで作ったのですが、 どうしても最後のほうがわかりません。 関数fとf2のreturn やメイン関数を変更すればいいとおもうのですが、どうやっても満足な結果が得られません。どなたか詳しい方教えてもらえない? #include<stdio.h> #include<math.h> double f(double t,double y1) { return y1; } double f2(double t2,double y2) { double w=5.0,k=0.5; return -w*w*f(t2,y2)-2.0*k*y2; } int main(void) { int count,bunkatu,npr; double t_0,T; double k1,k2,k3,k4; double y,t,dt,dt_2,dt_6; t_0 = 0.0; T=10.0; bunkatu=2000; npr=1000; y=1.0; printf("\n# [t_0,T]=[%6.4lf,%6.4lf] N=%d\n",t_0,T,bunkatu); t=t_0; dt=(T-t_0)/(double)bunkatu; dt_2=dt/2.0; dt_6=dt/6.0; for(count=0;count <=bunkatu;count++) { if(count %(bunkatu/npr)==0) printf("\n %5.2lf %10.5f",t,y); k1=f2(t,y); k2=f(t+dt_2,y+k1*dt_2); k3=f(t+dt_2,y+k2*dt_2); k4=f(t+dt,y+k3*dt); y+=(k1+2.0*k2+2.0*k3+k4)*dt_6; t+=dt; } printf("\n"); return 0; }

みんなの回答

  • Tacosan
  • ベストアンサー率23% (3656/15482)
回答No.5

1階微係数って記憶しなくてもいいんだっけ?

  • dra2jp
  • ベストアンサー率25% (18/72)
回答No.4

本当に教えてほしいんでしょうかね。 配布されたものらしきプログラムを投げつけ ルンゲクッタ法を用いたバネの単振動を解析する課題 としか説明せず みんながよくわからないなりに必死に読んでくれてるにもかかわらず、補足要求に たいしても 「初期値を決めて解析する方法」 としか言わず。 本当に教えてもらいたいなら プログラムの注釈をこまかくかいたり、 どこでどのような処理をするプログラムか、説明を書いたり、 物理の様々な式を計算しているなら 最低限計算している式を書くのが常識です。 私も読みましたがルンゲクッタ法がよくわからないのでプログラムをどう改良したらよいかもよくわかりませんでした。 詳しい補足要求お願いします。

  • mizuneko
  • ベストアンサー率16% (3/18)
回答No.3

物理わからないなりに読んでみましたが、 まず、fとf2の二種類あるのがよくわかりません。 また、"単振動を解析"だけでは何を求めるのかわかりませんが、 仮に求めるのが、時刻tにおける位置yのグラフであるなら、 関数f()は時刻tにおける速度sを返すはずです。 return y1; が速度でないということはわかりますが、どうでしょうか。 時刻tにおける速度sを求める計算式は、私にはわからないので ご自身で数学的に求めてください^^。

  • sirn
  • ベストアンサー率14% (8/55)
回答No.2

これはプログラミングするとき全般に言えることですが、 「他人に見せるコードにはわかりやすくコメントしましょう」を守ってください(^^; 私は大学に入って間もないのでルンゲクッタをよく知りません。 それで式と照らし合わそうと努力したのですが、大変な苦労です。 どの変数がルンゲクッタの式の何に対応しているのか教えてください。

  • sirn
  • ベストアンサー率14% (8/55)
回答No.1

まず、どういう処理をしているか書いてください。 何処が分からないのか分かりません。 ルンゲ=クッタ法の方法通りにやってください。

kodo_mo_dayo
質問者

補足

y''=-w^2y'-2.0ky で初期値を決めて解析する方法です。

関連するQ&A

  • 配列を用いたC言語プログラミングについて

    以下のルンゲクッタ法を用いたプログラムに配列などを使いさらに短くしたいのですが どのような方法が有りますか? #include <stdio.h> #include <math.h> double f1(double t1,double w,double x,double y,double z); double f2(double t1,double w,double x,double y,double z); double f3(double t1,double w,double x,double y,double z); double f4(double t1,double w,double x,double y,double z); //箱Aの関数 double g1(double t1,double a,double b,double c,double d); double g2(double t1,double a,double b,double c,double d); double g3(double t1,double a,double b,double c,double d); double g4(double t1,double a,double b,double c,double d); //箱Bの関数 int main(void) { double t1,w,x,y,z,a,b,c,d,dt,t1max,t2max,lam,gam,lat,dw,dx,dy,dz,da,db,dc,dd ; double k1[4],k2[4],k3[4],k4[4],l1[4],l2[4],l3[4],l4[4] ; ///宣言 t1 = 0.0; dt = 0.3; t1max = 40.0; //時間初期値 w = 200.0; x = 40.0; y = 30.0; z = 30.0; ///箱A初期値(w:感受性人口、x:潜伏人口、y:感染人口、z:隔離人口) a = 20.0; b = 8.0; c = 12.0; d = 10.0; ///箱B初期値(a:感受性人口、b:潜伏人口、c:感染人口,d:隔離人口) for(t1=0.0;t1<=t1max;t1+=dt) { k1[0]=dt*f1(t1,w,x,y,z); k1[1]=dt*f2(t1,w,x,y,z); k1[2]=dt*f3(t1,w,x,y,z); k1[3]=dt*f4(t1,w,x,y,z); k2[0]=dt*f1(t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0); k2[1]=dt*f2(t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0); k2[2]=dt*f3(t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0); k2[3]=dt*f4(t1+dt/2.0,w+k1[0]/2.0,x+k1[1]/2.0,y+k1[2]/2.0,z+k1[3]/2.0); k3[0]=dt*f1(t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0); k3[1]=dt*f2(t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0); k3[2]=dt*f3(t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0); k3[3]=dt*f4(t1+dt/2.0,w+k2[0]/2.0,x+k2[1]/2.0,y+k2[2]/2.0,z+k2[3]/2.0); k4[0]=dt*f1(t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]); k4[1]=dt*f2(t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]); k4[2]=dt*f3(t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]); k4[3]=dt*f4(t1+dt,w+k3[0],x+k3[1],y+k3[2],z+k3[3]); ///箱Aルンゲクッタ l1[0]=dt*g1(t1,a,b,c,d); l1[1]=dt*g2(t1,a,b,c,d); l1[2]=dt*g3(t1,a,b,c,d); l1[3]=dt*g4(t1,a,b,c,d); l2[0]=dt*g1(t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0); l2[1]=dt*g2(t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0); l2[2]=dt*g3(t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0); l2[3]=dt*g4(t1+dt/2.0,a+l1[0]/2.0,b+l1[1]/2.0,c+l1[2]/2.0,d+l1[3]/2.0); l3[0]=dt*g1(t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0); l3[1]=dt*g2(t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0); l3[2]=dt*g3(t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0); l3[3]=dt*g4(t1+dt/2.0,a+l2[0]/2.0,b+l2[1]/2.0,c+l2[2]/2.0,d+l2[3]/2.0); l4[0]=dt*g1(t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]); l4[1]=dt*g2(t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]); l4[2]=dt*g3(t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]); l4[3]=dt*g4(t1+dt,a+l3[0],b+l3[1],c+l3[2],d+l3[3]); ///箱Bルンゲクッタ w=w+((k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0); x=x+((k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0); y=y+((k1[2]+2.0*k2[2]+2.0*k3[2]+k4[2])/6.0); z=z+((k1[3]+2.0*k2[3]+2.0*k3[3]+k4[3])/6.0); a=a+((l1[0]+2.0*l2[0]+2.0*l3[0]+l4[0])/6.0); b=b+((l1[1]+2.0*l2[1]+2.0*l3[1]+l4[1])/6.0); c=c+((l1[2]+2.0*l2[2]+2.0*l3[2]+l4[2])/6.0); d=d+((l1[3]+2.0*l2[3]+2.0*l3[3]+l4[3])/6.0); } return 0; }

  • C言語のプログラムを作るのに困っています。

    C言語で、4次のルンゲクッタ法を用いて微分方程式を解くプログラムを作成したいのですが、始めたばかりの初心者で分からないところだらけなので教えてください。 わからないなりに、下のプログラムを作ってみましたがエラーになります。どうすればうまくいくのでしょうか。本当に始めたばかりなので、文?の組立もよくわかりません。できればわかりやすくお願いします。 このプログラムでは微分方程式dy/dx=x^3+4xy+y+2の解を求めようとしています(自分で適当に作った方程式なので、あるのかどうかわかりません)。 そのほかにもy'=2x^2+3xy-1や、y'=sinxcosy, y'=e^xなどいろいろな微分方程式でできるようにしたいのですが、どうやるのでしょうか。 #include <stdio.h> double f(double x,double y); int main() { double x,y,dx,xmax; double k1,k2,k3,k4; FILE *fp; if ((fp=fopen("result.txt","w"))==NULL) { printf("Cannot open result.txt \n"); return 0; } dx=0.1; //刻み幅 xmax=100.0; //繰り返し最大 x=1.0; //xの初期値 y=1.0; //yの初期値 f(x,y)=x^3+4xy+y+2;; //関数f(x,y)の定義 for(x=1.0;x<xmax;x+=0.1) { k1=dx*f(x,y); k2=dx*f(x+dx/2.0,y+k1/2.0); k3=dx*f(x+dx/2.0,y+k2/2.0); k4=dx*f(x+dx,y+k3); y=y+(k1+2.0*k2+2.0*k3+k4)/6.0; printf(" %f %lf \n",x,y); fprintf(fp," %f %lf \n",x ,y); } return 0; }

  • ルンゲクッタ法で数値解析(C言語)

    コンピュータ(プログラミング)のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。 単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。 書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました(C言語を使って)。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。 ---- #include <stdio.h> double f1(double t,double x,double v); double f2(double t,double x,double v); int main() { double t,x,v,dt,tmax; double k1[2],k2[2],k3[2],k4[2]; x=1.0; //位置の初期値 v=0.0; //速度の初期値 dt=0.01; //刻み幅 tmax=500.0; //繰り返し最大回数 FILE *output; output=fopen("output.dat","w"); for(t=0;t<tmax;t+=dt) { k1[0]=dt*f1(t,x,v); k1[1]=dt*f2(t,x,v); k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k3[0]=dt*f1(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k3[1]=dt*f2(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k4[0]=dt*f1(t+dt,x+k3[0],v+k3[1]); k4[1]=dt*f2(t+dt,x+k3[0],v+k3[1]); x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0; v=v+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0; fprintf(output,"%f %f %f\n",t,x,v); } fclose(output); return 0; } double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); } ---- このソースは単振動のもので、減衰振動のときは、最後の方の ---- double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); ---- の部分が変わるのだと思うのですが、よく分かりません。 減衰振動は、mx"=-kx-rx'で表されるので、この式を変形(?)したものが、入るのかな、という予想程度にしか分かりません。 ソースをどのように変えれば減衰振動を解析できるのでしょうか。 どなたか詳しい方、教えてください。お願いします。

  • C言語 二分法

    初投稿です。 お恥ずかしながらパソコンが苦手で、Cゲッが難しくてできません。 今回二分法です。 途中まではやったのですができません。 演習 0~1の乱数を12個発生させ,これらの平均をxi,yi とする。(s,tは乱数) xi=1/12(x1i+x2i+....x12i) yi=1/12(y1i+y2i+....y12i) このようなを1000個作り,(xi,yi)で散布図 を作りなさい。またx,yのそれぞれの平均を求 めよ。 この演習で #include<stdio.h> #include<stdlib.h> #include<math.h> #define eps 1.0e-5 float f(double x); void nibuin(void); int main () { int count; double a,b,m; count=0; printf("範囲の左の値を入力してください。\n"); scanf("%lf",&a); printf("範囲の右の値を入力してください。\n"); scanf("%lf",&a); if(count==1000){ printf("収束しませんでした。\n"); exit(1); } } while(!(fabs(a-b)<eps)); printf("解の値は%f\n収束するのに%d回かかりました。"m,count); } float f(double x) { reurn x*sin(x)+log(x); } まではできたのですが、 scanf("%lf",&a);とif(count==1000){の間に入る命令が打てません。 よろしくお願いします。

  • コンパイルした時に表示されません

    バネのつり合いから,オーバーシュートとその時の時間を求めるのですが,以下のようにプログラムを作り,オーバーシュートprintf("\n%lf\n\n",y1-y0);は出せましたが,一番下のprintf("t=\n%lf\n\n",(y1-y0)/(1+y1));の部分が画面に表示されません.わかる方ぜひ教えてください. #include<stdio.h> #define G 9.8 double m,K,H,A,N,t; double h=0.01; double f1(double y1,double y2){ return(y2); } double f2(double y1,double y2){ return(-K/m*y1+G-A*N/H/m*y2); } main(){ int i; double k[5][3]; double y0,y1,y2; y1=0; y2=0; printf("m,K,H,A,N="); scanf("%lf,%lf,%lf,%lf,%lf",&m,&K,&H,&A,&N); for(i=0;i<1000;i++){ k[1][1]=h*f1(y1,y2); k[1][2]=h*f2(y1,y2); k[2][1]=h*f1(y1+k[1][1]/2.0,y2+k[1][2]/2.0); k[2][2]=h*f2(y1+k[1][1]/2.0,y2+k[1][2]/2.0); k[3][1]=h*f1(y1+k[2][1]/2.0,y2+k[2][2]/2.0); k[3][2]=h*f2(y1+k[1][1]/2.0,y2+k[1][2]/2.0); k[4][1]=h*f1(y1+k[3][1],y2+k[3][2]); k[4][2]=h*f2(y1+k[3][1],y2+k[3][2]); if(k[1][1]+k[2][1]+k[3][1]+k[4][1]<0){ break; } y1=y1+k[1][1]/6.0+k[2][1]/3.0+k[3][1]/3.0+k[4][1]/6.0; y2=y2+k[1][2]/6.0+k[2][2]/3.0+k[3][2]/3.0+k[4][2]/6.0; } y0=m*G/K; printf("\n%lf\n\n",y1-y0); printf("t=\n%lf\n\n",(y1-y0)/(1+y1)); }

  • C言語課題(数値計算)お願いします。

    課題内容 ニュートン法を用いて√2の値を求められるプログラムを作る。 10ケタまで表示する。 書いてみたソースコード #include <stdio.h> #include <stdlib.h> #include <math.h> #define max 1000 //最大繰り返し回数 #define eps 1.0e-10 //収束条件 double f(double x); double df(double x); void newton(void); int main() { newton(); return 0; } void newton(void) { int count; double a,newa; count=0; printf("初期値を入力してください。\n"); scanf("%lf",&a); for(;;) { count++; newa=a-f(a)/df(a); if(fabs(newa-a)<eps) break; a=newa; if(count==max) { printf("収束しませんでした。\n"); exit(1); } } printf("解の値は %f\n収束するのに %d 回かかりました。\n", newa,count); } double f(double x) { return x*x-2; } double df(double x) { return 2x; } コンパイルできません。 間違いとどうすればよいかを教えていただけるとありがたいです!!!

  • (C言語)プログラムが動きません。

    下のC言語のプログラムを書いたのですが、動きません。上手く作動していれば画面に「結果をrunge1.csvに書き込みました」と表示されてテキストファイルが作成されるはずですが、プログラムを実行しても何のメッセージも表示されずファイルも作成されません。ですが、bcpadのエラーメッセージの欄には何も表示されず、どこを直せばいいのか分かりません。アドバイスを下さい。お願いします。 /*ルンゲ・クッタ法による微分方程式y'+y*y=xの解*/ #include<stdio.h> #include<math.h> #define N 100 int main(void) { char *file="runge1.csv"; double x0=0.0,x1=5.0,y0=0.0,x,y,h,k1,k2,k3,k4; int i; FILE *fp; h=(x0-x1)/N; x=x0; y=y0; fp=fopen(file,"w"); fprintf(fp,"x,Runge-y\n"); fprintf(fp,"%3.2f,%3.5f\n",x,y); for(i=1;i<=N;i++){ x+=h; k1=(x-y*y)*h; k2=(x-y*y+k1/2.0)*h; k3=(x-y*y+k2/2.0)*h; k4=(x-y*y+k3)*h; y=(x-y*y)+(k1+2*k2+2*k3+k4)/6.0; fprintf(fp,"%3.2f,%3.5f\n",x,y); } fclose(fp); printf("結果を%sに書き込みました。\n",file); return(0); }

  • コンパイルした時になぜ表示されないか

    バネのつり合いから,オーバーシュートとその時の時間を求めるのですが,以下のようにプログラムを作り,オーバーシュートは出せましたが,一番下のprintf("t=\n%lf\n\n",(y1-y0)/(1+y1));の部分が画面に表示されません.わかる方ぜひ教えてください. #include<stdio.h> #define G 9.8 double m,K,H,A,N,t; double h=0.01; double f1(double y1,double y2){ return(y2); } double f2(double y1,double y2){ return(-K/m*y1+G-A*N/H/m*y2); } main(){ int i; double k[5][3]; double y0,y1,y2; y1=0; y2=0; printf("m,K,H,A,N="); scanf("%lf,%lf,%lf,%lf,%lf",&m,&K,&H,&A,&N); for(i=0;i<1000;i++){ k[1][1]=h*f1(y1,y2); k[1][2]=h*f2(y1,y2); k[2][1]=h*f1(y1+k[1][1]/2.0,y2+k[1][2]/2.0); k[2][2]=h*f2(y1+k[1][1]/2.0,y2+k[1][2]/2.0); k[3][1]=h*f1(y1+k[2][1]/2.0,y2+k[2][2]/2.0); k[3][2]=h*f2(y1+k[1][1]/2.0,y2+k[1][2]/2.0); k[4][1]=h*f1(y1+k[3][1],y2+k[3][2]); k[4][2]=h*f2(y1+k[3][1],y2+k[3][2]); if(k[1][1]+k[2][1]+k[3][1]+k[4][1]<0){ break; } y1=y1+k[1][1]/6.0+k[2][1]/3.0+k[3][1]/3.0+k[4][1]/6.0; y2=y2+k[1][2]/6.0+k[2][2]/3.0+k[3][2]/3.0+k[4][2]/6.0; } y0=m*G/K; printf("\n%lf\n\n",y1-y0); printf("t=\n%lf\n\n",(y1-y0)/(1+y1)); }

  • 微分方程式のプログラムについて

    次の問題を解こうとしましたが、プログラムが動かなかったので実行できませんでした。今、非常に困っています。分かりにくいプログラムですみません。 問 温度Taの空気中に置かれている物体の温度Tは、次の常微分方程式に従う。    dT/dt=-k(T-Ta) t=0でT=100[℃]、Ta=20[℃]、k=1×10^-2[s^-1]として、物体の温度が60℃にまで冷却される時間を求めよ。厳密解とも比較すること。 自分が作ったプログラムは #include <stdio.h> #include <math.h> double k = 1e-2, tmpa=20.0; double f(double t, double tmp) { return -k*(tmp-tmpa); } int main(void) { double t0 = 0.0, tmp0 = 100, t1, tmp1, ttmp1; double k1, k2, k3, k4, h; printf("刻み幅="); scanf("%lf", h); t1 = t0+h; ttmp1=tmpa+(tmp0-tmpa)+exp(-k*(t0-t0)); printf("%8.4lf %8.4lf %8.4lf\n", t0, tmp0, ttmp1); while(t0 <= 100.0){ k1 = f(t0, tmp0); k2 = f(t0+0.5*h, tmp0+0.5*h*k1); k3 = f(t0+0.5*h, tmp0+0.5*h*k2); k4 = f(t0+h, tmp0+h*k3); tmp1 = tmp0+h*(k1+2.0*k2+2.0*k3+k4)/6.0; printf("%8.4lf %8.4lf %8.4lf\n", t1, tmp1, ttmp1); if(tmp1 <= 60.0) break; t0 = t1; tmp0 = tmp1; t1 = t0+h; } return 0; } あと厳密解はT=e^-kt(T0-Ta)+Ta求められるそうです。 変な間違いが多いと思いますが、このプログラムのどこが間違っているか教えてください。厳密解だけでもいいです。お願いします。

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

    正弦波(およびガウス性雑音)をフーリエ変換(離散)→逆フーリエ変換するというプログラムを組みました。正弦波をフーリエ変換すると実部は2回ピークがくるはずなのですが、すべて「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; }

専門家に質問してみよう