Q&A
前提・実現したいこと
c言語でルンゲクッタ法を用いて出したデータをプロットしたいです。
プログラム内でファイルに書き出す方法で行っています。
発生している問題・エラーメッセージ
gnuplotを用いてプロットしたいのですが、横軸 時間、縦軸 摩擦力 としたいです。
その場合、FILE *fp;...の部分の書き方とそのコードの位置を教えていただきたいです。
最終的には、θなどを変えて摩擦力を見ていく予定です。
下記のプログラムでプロットした時のgnuplotでのWariningメッセージです。
gnuplot> plot "run.dat" w p pt 7 title "friction " Warning: empty x range [1.001:1.001], adjusting to [0.99099:1.01101] Warning: empty y range [6:6], adjusting to [5.94:6.06] Warning: slow font initializationqt.qpa.fonts: Populating font family aliases took 1657 ms. Replace uses of missing font family "Sans" with one that exists to avoid this cost.
該当のソースコード
#include<stdio.h> #include<math.h> double f1(double t, double x1, double x2, double x3, double vold, double w, double g); double f2(double t, double x1, double x2, double x3, double vold, double w, double g, int i); int main(){ int i,n; //double alpha; double dt=0.001,tmax,k=23,t=0.0,fk,fs,theta,fz,a,b,g,c; int Nsd[10]; double kx1[10],kv1[10],kx2[10],kv2[10],kx3[10],kv3[10],kx4[10],kv4[10],kw1[10],kw2[10],kw3[10],kw4[10],kg1[10],kg2[10],kg3[10],kg4[10]; double xold[10],xnew[10],vold[10],vnew[10],wold[10],wnew[10],gold[10],gnew[10]; double f[10]; double alpha[10]; n=10; //ブロックの個数 fs=0.7; //静摩擦係数 fk=0.45; //動摩擦係数 fz=60.0; //全荷重 tmax=1.001; theta=0.0; //初期条件 xold[0]=0.0; //xの初期値=0 for(i=0;i<=n-1;i+=1) { xold[i+1]=xold[i]+1.0; } //wold[0]=1.0; for(i=0;i<=n-1;i+=1) { wold[i]=0.0;//全ブロックの荷重の初期値 } for(i=0;i<=n-1;i+=1) { gold[i]=0; //全部のブロックに横向きの力 } for(t=0.0;t<=tmax;t+=dt) { for(i=0;i<=n-1;i++){ Nsd[i]=0; //静的 } for(i=0;i<=n-1;i++){ if(i==0){ if(fs*wold[i] < gold[i]+k*(xold[i+1]-xold[i])) { Nsd[i]=1; //動的 f[i]=fk*wold[i]; //printf("%lf\n",f[i]); } } else if(i==n-1){ if(fs*wold[i] < gold[i]+k*(-xold[i]+xold[i-1])) { Nsd[i]=1; f[i]=fk*wold[i]; //printf("%d %lf\n",i,f[i]); } } else{ if(fs*wold[i] < gold[i]+k*(xold[i+1]-2*xold[i]+xold[i-1])) { Nsd[i]=1; f[i]=fk*wold[i]; //printf("%d %lf\n",i,f[i]); } } } /* for(i=0;i<=n-1;i++) { if(Nsd[i]==0) { f[i]=gold[i]; } } */ for(i=0;i<=n-1;i++){ kx1[i]=dt*f1(t,xold[i-1],xold[i],xold[i+1],vold[i],wold[i],gold[i]); kv1[i]=dt*f2(t,xold[i-1],xold[i],xold[i+1],vold[i],wold[i],gold[i],i); kw1[i]=0; kg1[i]=0; if(Nsd[i]==0){ kx1[i]=0; kv1[i]=0; } } for(i=0;i<=n-1;i++){ kx2[i]=dt*f1(t+dt*0.5,xold[i-1]+kx1[i-1]*0.5,xold[i]+kx1[i]*0.5,xold[i+1]+kx1[i+1]*0.5,vold[i]+kv1[i]*0.5,wold[i]+kw1[i]*0.5,gold[i]+kg1[i]*0.5); kv2[i]=dt*f2(t+dt*0.5,xold[i-1]+kx1[i-1]*0.5,xold[i]+kx1[i]*0.5,xold[i+1]+kx1[i+1]*0.5,vold[i]+kv1[i]*0.5,wold[i]+kw1[i]*0.5,gold[i]+kg1[i]*0.5,i); kw2[i]=0; kg2[i]=0; if(Nsd[i]==0){ kx2[i]=0; kv2[i]=0; } } for(i=0;i<=n-1;i++){ kx3[i]=dt*f1(t+dt*0.5,xold[i-1]+kx2[i-1]*0.5,xold[i]+kx2[i]*0.5,xold[i+1]+kx2[i+1]*0.5,vold[i]+kv2[i]*0.5,wold[i]+kw2[i]*0.5,gold[i]+kg2[i]*0.5); kv3[i]=dt*f2(t+dt*0.5,xold[i-1]+kx2[i-1]*0.5,xold[i]+kx2[i]*0.5,xold[i+1]+kx2[i+1]*0.5,vold[i]+kv2[i]*0.5,wold[i]+kw2[i]*0.5,gold[i]+kg2[i]*0.5,i); kw3[i]=0; kg3[i]=0; if(Nsd[i]==0){ kx3[i]=0; kv3[i]=0; } if(i==n-1){ kx3[i]=dt*f1(t+dt*0.5,xold[i-1]+kx2[i-1]*0.5,xold[i]+kx2[i]*0.5,0,vold[i]+kv2[i]*0.5,wold[i]+kw2[i]*0.5,gold[i]+kg2[i]*0.5); kv3[i]=dt*f2(t+dt*0.5,xold[i-1]+kx2[i-1]*0.5,xold[i]+kx2[i]*0.5,0,vold[i]+kv2[i]*0.5,wold[i]+kw2[i]*0.5,gold[i]+kg2[i]*0.5,i); kw3[i]=0; kg3[i]=0; if(Nsd[i]==0){ kx3[i]=0; kv3[i]=0; } } } for(i=0;i<=n-1;i++){ kx4[i]=dt*f1(t+dt,xold[i-1]+kx3[i-1],xold[i]+kx3[i],xold[i+1]+kx3[i+1],vold[i]+kv3[i],wold[i]+kw3[i],gold[i]+kg3[i]); kv4[i]=dt*f2(t+dt,xold[i-1]+kx3[i-1],xold[i]+kx3[i],xold[i+1]+kx3[i+1],vold[i]+kv3[i],wold[i]+kw3[i],gold[i]+kg3[i],i); kw4[i]=0; kg4[i]=0; if(Nsd[i]==0){ kx4[i]=0; kv4[i]=0; } } for(i=0;i<=n-1;i++){ xnew[i]=xold[i]+(kx1[i]+2.0*kx2[i]+2.0*kx3[i]+kx4[i])/6.0; vnew[i]=vold[i]+(kv1[i]+2.0*kv2[i]+2.0*kv3[i]+kv4[i])/6.0; xold[i]=xnew[i]; vold[i]=vnew[i]; } alpha[9]=10.0; //alpha[0]=5.5; for(i=0;i<=n-1;i++){ //alpha[i+1]=alpha[i]+1.0; //gnew[i]=alpha[i]*t; gnew[i]=alpha[i]*t; gold[i]=gnew[i]; //printf("%lf",gold[i]); if(Nsd[i]==0){ f[i]=gold[i]; } } for(i=0;i<=n-1;i+=1) { //theta+=0.001; wold[i]=(fz/n)*(1+(2*(i+1)-n-1)*theta/(n-1)); //荷重 } b=k*(xold[2]-2*xold[1]+xold[0]); c=(f[0]+f[1]+f[2]+f[3]+f[4]+f[5]+f[6]+f[7]+f[8]+f[9])/10;//摩擦力 g=(gold[0]+gold[1]+gold[2]+gold[3]+gold[4]+gold[5]+gold[6]+gold[7]+gold[8]+gold[9])/10;//駆動力 //printf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",t,f[0],f[1],f[2],f[3],f[4],f[5],f[6],f[7],f[8],f[9]); //printf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",t,xold[0],xold[1],xold[2],xold[3],xold[4],xold[5],xold[6],xold[7],xold[8],xold[9]); //printf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",t,wold[0],wold[1],wold[2],wold[3],wold[4],wold[5],wold[6],wold[7],wold[8],wold[9]); //printf("%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",t,gold[0],gold[1],gold[2],gold[3],gold[4],gold[5],gold[6],gold[7],gold[8],gold[9]); //printf("%lf %d %d %d %d %d %d %d %d %d %d\n",t,Nsd[0],Nsd[1],Nsd[2],Nsd[3],Nsd[4],Nsd[5],Nsd[6],Nsd[7],Nsd[8],Nsd[9]); //printf("%lf %lf\n",t,theta); //printf("%lf %lf\n",t,g); printf("%lf %lf\n",t,c); //printf("%lf %lf\n",t,f[8]); //printf("%lf %lf\n",t,b); } FILE *fp; fp=fopen("run.dat","w"); fprintf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",t,wold[0],wold[1],wold[2],wold[3],wold[4],wold[5],wold[6],wold[7],wold[8],wold[9]); } double f1(double t,double x1, double x2, double x3, double vold, double w, double g) { return (vold); } double f2(double t,double x1, double x2, double x3, double vold, double w, double g, int i) { double fk=0.45,fs=0.7; double a2; int n=10; if(vold==0){ a2=(x1-2.*x2+x3)+g; //a2=0.; if(i==0){ a2=(x3-x2)+g; //a2=0.; } if(i==9){ a2=(x1-x2)+g; //a2=0.; } } else if(vold>0){ a2=(x1-2.*x2+x3)-fk*w+g; // printf("%lf",a2); if(i==0){ a2=(x3-x2)-fk*w+g; } if(i==9){ a2=(x1-x2)-fk*w+g; // printf("%lf",a2); } } else{ a2=(x1-2.*x2+x3)+fk*w+g; //a2=0.; if(i==0){ a2=(x3-x2)+fk*w+g; //a2=0.; } if(i==9){ a2=(x1-x2)+fk*w+g; //a2=0.; } } return (a2); }
上記では荷重としているwについて書き込もうとしていますが無視してください。
補足情報（FW/ツールのバージョンなど）
OS:Big Sur
言語：C (Xcode)
プロットはgnuplotを使っています。
10個のブロックをバネで繋げたモデルで、走る際の荷重の掛け方による摩擦の変化を見るプログラムです。
つま先、足裏全体、踵の３パターンをθなどの値を変えて調べたいと思っています。
知識不足ですが、お願い致します。
回答1件
0
ベストアンサー
現状では以下の様になっています。
c
1 for(t=0.0;t<=tmax;t+=dt) 2 { 3 : 4 printf("%lf %lf\n",t,c); 5 : 6 } 7 FILE *fp; 8 fp=fopen("run.dat","w"); 9 fprintf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",t,wold[0],wold[1],wold[2],wold[3],wold[4],wold[5],wold[6],wold[7],wold[8],wold[9]);
既にループ内で printf("%lf %lf\n",t,c); としているので、これを fprintf() に変更します。
c
1 FILE *fp = fopen("run.dat", "w"); 2 for(t=0.0;t<=tmax;t+=dt) 3 { 4 : 5 // printf("%lf %lf\n",t,c); 6 fprintf(fp, "%lf %lf\n", t, c); 7 : 8 } 9 fclose(fp); 10 11 return 0; 12}
投稿2021/11/22 08:12
総合スコア19865
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2021/11/23 00:20
2021/11/23 00:24
2021/11/23 00:55
2021/11/23 01:11 編集
2021/11/23 01:22
2021/11/23 01:37
2021/11/23 05:18
2021/11/23 05:23
2021/11/23 05:26