• 締切済み

数値解析法

このHeun法のプログラムをRunge-Kutta法にするにはどうしたらいいですか？ #include <stdio.h> #include <math.h> double f1(double y) { return y; } double f2(double y) { return -4*y; } int main(){ double a=0; double b; int m=10; int n; double h; double x,y; int k; double e; double f; double k1,k2; printf("Heun法計算例：y=e^x, y=1/e^4x\n\n"); // y = e^x b = 1; for(n=100;n<=10000;n*=100){ h = (b-a)/n; printf("y' = y: h(=dx) = %.1e (y=e^x)\n",h); x = a; y = 1; for(k=0;k<=n;k++) { x = k*h; if(k%(n/m)==0) { f = exp(x); e = fabs(y-f); printf("x=%.2f, y=%f, e^x=%f er=%.0e\n",x,y,f,e); } // Heun's method k1 = h*f1(y); k2 = h*f1(y+k1); y += (k1+k2)/2; } }

みんなの回答

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

どうしても困ったら名前だけ変えて「Runge-Kutta」っていいはればいいのでは? あなたに十分なバックグラウンドがあれば通せるかもしれませんよ.

  1. カテゴリ
  2. [技術者向] コンピューター
  3. プログラミング・開発
  4. C・C++・C#

関連するQ&A

  • C言語のプログラム

    C言語で2つの微分方程式をEuler法、Heun法、Runge-Kutta法により求めるプログラムを作りたい。ただし、初期条件はx=0,y=1とする。また、間隔Δxを変えたときの解の変化を調べたい。 Euler法のプログラムはどうにか分かったのですが、Heun法、Runge-Kutta法のプログラムがわかりません。 Euler法のプログラム #include_<stdio.h> #include_<math.h> int_main(){ __double__a=0; __double__b; __int_____m=10; __int_____n; __double__h; __double__x,y; __double__dydx; __int_____k; __double__e; __double__f; __printf("オイラー法計算例：y=e^x,_y=1/e^4x\n\n"); __//_y_=_e^x __b_=_1; __for(n=100;n<=10000;n*=100){ ____h_=_(b-a)/n; ____printf("y'_=_y:_h(=dx)_=_%.1e_(y=e^x)\n",h); ____x_=_a;_y_=_1; ____for(k=0;k<=n;k++)_{ ______x_=_k*h; ______if(k%(n/m)==0)_{ ________f_=_exp(x); ________e_=_fabs(y-f); ________printf("x=%.2f,_y=%f,_e^x=%f_er=%.0e\n", ________x,y,f,e); ______} ______dydx_=_y; ______y_=_y_+_dydx*h; ____} __} __printf("\n"); __//_y_=_1/e^4x __b_=_4; __for(n=100;n<=10000;n*=100){ ____h_=_(b-a)/n; ____printf("y'_=_-4y:_h(=dx)_=_%.1e_(y=1/e^4x)\n",h); ____x_=_a;_y_=_1; ____for(k=0;k<=n;k++)_{ ______x_=_k*h; ______if(k%(n/m)==0)_{ ________f_=_exp(-4*x); ________e_=_fabs(y-f); ________printf("x=%.2f,_y=%f,_1/e^4x=%f_er=%.0e\n", ________x,y,f,e); ______} ______dydx_=_-4*y; ______y_=_y_+_dydx*h; ____} __} __return_0; } 分かる方がいましたら、回答よろしくお願いします。

  • ルンゲクッタのプロクラム実行できません

    プログラムの実行ができません。 どこを直したらよいのか教えて下さい。 よろしくお願いします。 #include <stdio.h> #include <stdlib.h> #include <math.h> double *dvector(long i,long j); void free_dvector(double *a,long i); double func(double x,double y); double *rk4(double y0,double a,double b, int n,double(*f)()); int main(void) { double *y,h,a=0.0, b=1.0, y0=1.0; int i,n; printf("分割数を入力して下さい------->"); scanf("%d",&n); y=dvector(0,n); y=rk4(y0,a,b,n,func); h=(b-a)/n; for(i=0;i<=n;i++) { printf("x=%f \t y=%f \n",a+i*h,y[i]); } return 0; } double *rk4(double y0,double a, double b,int n,double (*f)()) { double k1,k2,k3,k4,h,*y,x; int i; y=dvector(0,n); h=(b-a)/n; y[0]=y0; x=a; for(i=0;i<n;i++) { k1=f(x,y[i]); k2=f(x+h/2.0,y[i]+h*k1/2.0); k3=f(x+h/2.0,y[i]+h*k2/2.0); k4=f(x+h,y[i]+h*k3); y[i+1]=y[i]+h/6.0*(k1+2.0*k2+2.0*k3+k4); x +=h; } return y; free_dvector(y,0); }

  • Xcodeでのc言語のプログクムについて。

    Xcodeでのc言語のプログクムについて。 ルンゲークッタ法のプログラムを書いているのですがビルドエラーになってしまい、原因が分かりません。 　ビルド結果は、 "_func",referenced from: _main in main.o "_dvector",referenced from: _main in main.o _rk4 in main.o symbol(s) not found collect2: ld returned 1 exit status 　プログラムは、 #include <stdio.h> #include <stdlib.h> double *dvector(long i, long j); /* ベクトル量領域の確保*/ void free_dvector(double *a, long i); /*　領域の解放*/ double func(double x, double y); /* 関数の定義 */ /* ルンゲークッタ法 */ double *rk4( double y0, double a, double b, int n, double (*f)()); int main (void) { double *y, h, a=0.0, b=1.0, y0=1.0 ; int i, n; printf("分割式を入力してください-->"); scanf("%d",&n); y = dvector( 0, n ); /*領域の確保*/ y = rk4( y0, a, b, n, func ); /*ルンゲクッタ法*/ /*　結果の表示 */ h = (b-a)/n ; /* 刻み幅　*/ for (i = 0 ; i <= n ; i++) { printf("x=%f \t y=%f \n", a+i*h, y[i] ); } return 0; } /* ルンゲクッタ法　*/ double *rk4( double y0, double a, double b, int n, double (*f)() ) { double k1, k2, k3, k4, h, *y, x; int i; y = dvector( 0, n ); /* y[0,1,....,n]の確保　*/ h = (b-a)/n ; /* 初期値の設定　*/ y[0] = y0; x = a; /* ルンゲクッタ法　*/ for ( i = 0 ; i < n ; i++) { k1 = f(x,y[i]); k2 = f(x+h/2.0, y[i]+h*k1/2.0); k3 = f(x+h/2.0, y[i]+h*k2/2.0); k4 = f(x+h, y[i]+h*k3); y[i+1] = y[i] + h/6.0 * ( k1 + 2.0*k2 + 2.0*k3 + k4 ); x += y[i]; } return y; free_dvector(y, 0); /* 領域の解放　*/ }

  • 数値解析の矩形法について

    区間分割数nを10から100まで10ずつ増やして計算値を求め、面積の計算誤差が区間分割数によってどのように変化するかを求める矩形法プログラムを下記に作成したのですが、 計算誤差を求める式で ・計算誤差＝（計算値ー真値）/真値×100 （％） （真値は0.68269） と計算するのですが、プログラムでは分母の真値×100 （％）を（）を付ける場合の計算の答えと（）を付けない計算の答えとが全然違います。どうしてこのようなことが起こるのですか？ また、この計算誤差の求め方は（）を付ける場合とつけない場合のどちらが正しいのですか？ public class Kukei { static double f(double x) { // ここに任意の被積分関数を記述 double y = Math.exp(- x * x / 2) / Math.sqrt(2.0 * Math.PI); return y; } public static void main(String[] args) { double a = - 1.0, b = 1.0; // 積分範囲 int n = 10; // 区間分割数 double suti= ((n+10)-0.68269)/0.68269*100;　//真値 for(int k=1; k <= 10; k++){ double h = (b - a) / (double)n; // 分割幅 double s = 0.0; n=k*10; for (int i=0; i < n; i++) { s += f(a + i * h); } s *= h; suti= (s-0.68269)/0.68269*100; System.out.println("区間分割数　＝" + n); System.out.println("矩形法による計算値　＝" + s); System.out.println("矩形法による計算誤差　＝" +suti+"\n"); } } }

  • 二分法のプログラムについて

    下の用なプログラムを作ったのですがどうしても正しい答えを導くことができません。自分でもいろいろ調べてみましたがわかりません。誰かご教授宜しくお願いします。 #include<stdio.h> #include<stdlib.h> #define MAX 10 int n , count; double c[MAX+1]; double a,b,e; void nyuuryoku(void) { int i; printf("nの入力>"); scanf("%d",&n); if(n>MAX){printf("最大次数を超えている");exit(1);} else if(n<0){printf("nが負");exit(2);} else{for(i=0;i<=n;i++){printf("係数の値>");scanf("%lf",&c[i]);} }} double f(double x) {double y; int i; y = c[0]; for(i=1;i<=n;i++){ y=y*x+c[i];} return y; } void hani(void){ printf("aの値>");scanf("%lf",&a); printf("bの値>");scanf("%lf",&b); printf("eの値>");scanf("%lf",&e); if(e<=0){printf("eが0または負"); exit(3);} if(f(a)==0){printf("%f",f(a)); exit(4);} if(f(b)==0){printf("%f",f(b)); exit(5);} if(f(a)*f(b)>0){printf("初期値異常"); exit(6);}} double nibun(void) {double c; if(b>a){ while(b-a>e){ count++; c=(a+b)/2; if(f(c)==0){ return c;} if(f(a)*f(c)<0){b=c;} if(f(b)*f(c)<0){a=c;} } return a;} if(a>b){ while(a-b>e){ count++; c=(a+b)/2; if(f(c)==0){ return c;} if(f(b)*f(c)<0){a=c;} if(f(a)*f(c)<0){b=c;} } return a;} } void syutsuryoku(double x){ printf("x=%lf\n",x); printf("f(x)=%lf\n",f(x)); printf("繰り返し回数=%d\n",count); } int main(void){ double ans; count=0; nyuuryoku(); hani(); ans = nibun(); syutsuryoku(ans); }

  • ルンゲクッタ法で数値解析（Ｃ言語）

    コンピュータ（プログラミング）のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。 単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。 書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました（Ｃ言語を使って）。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。 ---- #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'で表されるので、この式を変形（？）したものが、入るのかな、という予想程度にしか分かりません。 ソースをどのように変えれば減衰振動を解析できるのでしょうか。 どなたか詳しい方、教えてください。お願いします。

  • 二分法のプログラム

    関数x^３－７x^２－６x＋２を二分法で解くプログラムを作ったのですが、エラーが出てコンパイルできません。訂正箇所を教えて下さい。 宜しくお願い致します。 #include<stdio.h> #include<math.h> #define EPSILON 0.1E-5 #define TURE 1 #define FALSE 0 int kansu(int x); void Nibunho(left,right,sol,flag); double left,right; int flag; int main(void) { printf("区間の左端と右端は？\n"); scanf("%lf %lf",&left,&right); flag=FALSE; Nibunho(left,right,&root,&flag); if(flag) printf("解 = %e (繰り返し回数 = %d)\n",root,k); else { printf("入力した範囲で解は求まりませんでした。\n"); printf("f(%e) = %e \n",root,k); } return 0; } int kansu(int x) { int f(double x) f(x)=x*x*x-7.0*x*x-6.0*x+2.0; return(f(x)); } void Nibunho(left,right,sol,flag) { double left,right,*sol; int *flag; double a,b,c,fa,fb,fc; k=0; a=left; b=right; do { k++; c=(a+b)/2.0; fc=f(c); fa=f(a); fb=f(b); if(fabs(fc)<1.0e-10) { a=c; b=c; *flag=TRUE; } else { if( (fa * fc < 0.0) || (fb * fc < 0.0) ) { *flag = TRUE; if( (fa*fc) < 0.0 ) b=c; else a=c; } else { if( fabs(fa) > fabs(fb) ) a=c; else b=c; } } } while((b-a)>EPSILON) *sol=(a+b)/2.0; }

  • ２分法で方程式の複数の解を自動的に求めるには？

    ２分法を授業で習い、アルゴリズムは理解できました。 そして、２分法で方程式の解を求めるプログラムも完成したのですが、 ３次方程式などの全ての解を２分法で求める場合、本来ならば、 グラフなどを描き、適切な初期値を自分で変更していく必要があると思います。 しかし、方程式の複数の解を自動的に求めたいのです。 もし、２分法で方程式の複数の解を自動的に求めるアルゴリズムがあれば教えていただけないでしょうか。 よろしくお願いします。 ※ 私が作成したプログラム（C言語）も載せておきます。 #include <stdio.h> #include <math.h> double f(double x) { return ((x-1.0)*(x-2.0)*(x-3.0)); } int main () { int i; double a,b,x; double gosa = 1.0e-14; printf("aの値を入力してください。\n"); scanf("%lf",&a); printf("bの値を入力してください。\n"); scanf("%lf",&b); if(f(a)*f(b) >0) printf("aとbの範囲の中に適切な解が存在しません。\n"); while(1){ x = (a+b)/2.0; printf("%.14f\n",x); if(fabs(b-a)<gosa) break; if(f(a)*f(x)<=0){ b = x; } else{ a = x; } } return(0); }

  • 大学の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; }

  • モンテカルロ法を用いた実験

    モンテカルロ法を用いてy=x^2とx軸で囲まれた部分の面積の近似値を求める．これを利用して平均μ，標準偏差σ，を求めたところ下記の写真のようになった．点数nに対して標準偏差は何次であるかを，中心極限定理から求められる次数，実験結果から導かれる次数のそれぞれで求めなさい．次数とはどのように求められるのでしょうか？使用したプログラムを下に書いておきます． #include <stdio.h> #include <time.h> #include <stdlib.h> #include <math.h> double f(double x); int main(void) { int M = 1000; int i,j,n,c,k; int tmp[6] = { 20, 50, 100, 200, 500, 1000 }; double x,y,Y, sumY, sumY2, muY, muY2, sigmaY; srand((unsigned int)time(NULL)); for (k = 0; k < 6; k++){ sumY = 0; sumY2 = 0; for (j = 0; j < M; j++){ c = 0; n = tmp[k]; for (i = 0; i < n; i++){ x = (1.0* (double)rand()) / ((double)RAND_MAX + 1.0); y = (1.0* (double)rand()) / ((double)RAND_MAX + 1.0); if (y <= f(x)){ c++; } } Y = ((double)c) / (double)n; sumY = sumY + Y; sumY2 += Y*Y; } muY = sumY / (double)M; muY2 =sumY2 / (double)M; sigmaY = sqrt(muY2 -(muY*muY)); printf("%d %lf %lf\n", n, muY, sigmaY); } return 0; } double f(double x) { double r = x * x; return r; }

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

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

質問する