Q&A
時刻t と そのときの摩擦 の組をテキストファイルに書き出せば、gnuplot で描画できるんですよね?
前提・実現したいこと
c言語でデータをプロットしたいです。
プログラム内でファイルに書き出す方法で行っています。
発生している問題・エラーメッセージ
真ん中あたりのコメントのprintfを用いて、横軸時間、縦軸摩擦力のグラフをプロットしたいです。
fk,fsは摩擦係数、woldは荷重でf[i]=fk*woldで摩擦力としています。
cは各ブロックの摩擦力の平均をとったものです。
最終的には、θなどを変えて摩擦力を見ていく予定です。
運動方程式は以下です。
Xが位置、Fが摩擦力、Gが駆動力です。
荷重の式は以下です。
これが駆動力です。
該当のソースコード
c #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,phi,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; phi=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; //全部のブロックに横向きの力 } FILE *fp = fopen("run.dat", "w"); 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]/n)*(1+(2*(i+1)-n-1)*phi/(n-1)*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]); //fprintf(fp,"%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",i,wold[0],wold[1],wold[2],wold[3],wold[4],wold[5],wold[6],wold[7],wold[8],wold[9]); //okfprintf(fp,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",theta,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]); //fprintf(fp,"%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); //fprintf(fp,"%d %lf\n",i,gnew[i]); **fprintf(fp,"%lf %lf\n",t,c);** //fprintf(fp,"%lf %lf\n",t,f[9]); //printf("%lf %lf\n",t,b); } fclose(fp); return 0; } 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); }
上のグラフを再現したいです。θ=0の摩擦力（平均）の時間変化です。
お願いいたします。
補足情報（FW/ツールのバージョンなど）
OS:Big Sur
言語：C (Xcode)
プロットはgnuplotを使っています。
10個のブロックをバネで繋げたモデルで、走る際の荷重の掛け方による摩擦の変化を見るプログラムです。
つま先、足裏全体、踵の３パターンをθなどの値を変えて調べたいと思っています。
知識不足ですが、お願い致します。
コードはその纏まりの前後を ``` の行で挟めば # で始まる行もちゃんとコードとして表示されます。
その際「ここに言語名を入力」という文字列を「c」に置き換えて頂ければ、コードを c 言語のソースとして扱ってくれます。
epistemeさん プロットはできます。画像のようなグラフにしたいのですが、どこを直すべきでしょうか。
jimbeさん ソースコードを貼る際に誤って消してしまっていたようです。ありがとうございます。
ご質問は編集できますので、編集して頂けると助かります。
ご質問自体につきまして、「C 言語で gnuplot によるプロットが出来ない」ということと、「”10個のブロックをバネで繋げたモデルで、走る際の荷重の掛け方による摩擦の変化を見る” プロットが出来ない」ということは、全く異なる分野です。
前者は C 言語の知識で分かりますが、後者は C 言語の知識では分かりません。
「式ははっきり分かっているが C 言語でその式をどう書くのか分からない」ということであれば、数学の知識と C 言語の知識で分かるのかもしれませんので、その式をご提示頂く必要があると思います。
回答2件
0
ベストアンサー
思うのですが、荷重との比較には絶対値を使うのではないでしょうか？
c
1for(i=0;i<=n-1;i++){ 2 if(i==0){ 3 //if(fs*wold[i] < gold[i]+k*(xold[i+1]-xold[i])) 4 if(fs*wold[i] < gold[i]+k*fabs(xold[i+1]-xold[i])) 5 { 6 Nsd[i]=1; //動的 7 f[i]=fk*wold[i]; 8 } 9 } 10 else if(i==n-1){ 11 // if(fs*wold[i] < gold[i]+k*(-xold[i]+xold[i-1])) 12 if(fs*wold[i] < gold[i]+k*fabs(-xold[i]+xold[i-1])) 13 { 14 Nsd[i]=1; 15 f[i]=fk*wold[i]; 16 } 17 } 18 else{ 19 // if(fs*wold[i] < gold[i]+k*(xold[i+1]-2*xold[i]+xold[i-1])) 20 if(fs*wold[i] < gold[i]+k*(fabs(xold[i+1]-xold[i])+fabs(-xold[i]+xold[i-1]))) 21 { 22 Nsd[i]=1; 23 f[i]=fk*wold[i]; 24 } 25 } 26}
また、荷重の値の更新部分は以下の様になるのではないでしょうか。
c
1 for(i=0;i<=n-1;i+=1){ 2 theta+=0.001; 3 //wold[i]=(fz/n)*(1+(2*(i+1)-n-1)*theta/(n-1)); //荷重 4 wold[i]=(fz/n)*(1+(2*(i+1)-(n-1))*theta/(n-1)); 5 }
2021/12/12 記載
配列に対する範囲外アクセスが発生しています。以下は修正したコードとの差分です。
diff
1@@ -25,7 +25,7 @@ 2 3 //初期条件 4 xold[0]=0.0; //xの初期値=0 5-for(i=0;i<=n-1;i+=1) 6+for(i=0;i<n-1;i+=1) 7 { 8 xold[i+1]=xold[i]+1.0; 9 } 10@@ -79,7 +79,7 @@ 11 } 12 } */ 13 14-for(i=0;i<=n-1;i++){ 15+for(i=1;i<n-1;i++){ 16 kx1[i]=dt*f1(t,xold[i-1],xold[i],xold[i+1],vold[i],wold[i],gold[i]); 17 kv1[i]=dt*f2(t,xold[i-1],xold[i],xold[i+1],vold[i],wold[i],gold[i],i); 18 kw1[i]=0; 19@@ -90,7 +90,7 @@ 20 } 21 } 22 23-for(i=0;i<=n-1;i++){ 24+for(i=1;i<n-1;i++){ 25 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); 26 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); 27 kw2[i]=0; 28@@ -100,7 +100,7 @@ 29 kv2[i]=0; 30 } 31 } 32-for(i=0;i<=n-1;i++){ 33+for(i=1;i<n-1;i++){ 34 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); 35 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); 36 kw3[i]=0; 37@@ -120,7 +120,7 @@ 38 } 39 } 40 } 41-for(i=0;i<=n-1;i++){ 42+for(i=1;i<n-1;i++){ 43 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]); 44 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); 45 kw4[i]=0;
しかし、この修正を行うと kx1[0], kx1[n-1], kx2[0], kx2[n-1], kx3[0], kx3[n-1], kx4[0], kx4[n-1] の値が未設定のままになりますので適宜埋める必要があります。
投稿2021/12/11 16:43
編集2021/12/13 02:52
総合スコア19872
回答ありがとうございます。修正後、コンパイル自体はできてプロットもできました。
未設定の値をどう決めたらいいかわかりません。
そこが最も重要なのだと思います。しかし、kx1[0] や kx1[n-1] などの値をどの様に計算するのが適切なのかは b6l9kk さんにしか解らないことなのです。。。
そうですよね、がんばります。配列の範囲外アクセスはエラーメッセージなどでわかるのでしょうか？
gcc の memory sanitizer という機能を使います。ただ、提示されているソースコードを先程読み返していたところ、回答した修正を行わなくても問題の無い事が判りました。f1 関数は vold を返すだけですし、f2 関数内では適切に計算が行われています。「適切に」という意味は、f2 関数内で if (i==0) とか if (i==9) としていて、実際には配列の範囲外の値を使ってはいないからです。ただ、f2 関数を呼び出す際に xold[i-1] や xold[i+1] としている部分では範囲外を参照していますので、そこは直しておく方が良いかと思います。
わかりました。ありがとうございます。
なぜ絶対値にするのでしょうか。
「荷重が負の値になることはないのでは？」と思うからです（この推測が正しいのかどうかは私には判断がつきません）。それと荷重で気が付いたのですが、荷重の値の更新部分に間違いがある様に思われます。回答に書きましたのでご確認下さい。
おっしゃる通りです。見落としていました。
参考にさせていただきます。ありがとうございます！
0
質問文にあるプログラムで出力されたファイルが gnuplot で表示されることは確認しました。
c言語：データのプロットがうまくいかない
例示されたグラフのイメージと同じにならない、という質問でしょうか？
10個のブロックをバネで繋げたモデルで、走る際の荷重の掛け方による摩擦の変化を見るプログラムです。
つま先、足裏全体、踵の３パターンをθなどの値を変えて調べたいと思っています。
アルゴリズムの話ということでしたら、そのあたりの知識は無いのでお答えしかねます。。
投稿2021/12/11 09:53
編集2021/12/11 09:57
総合スコア4633
あなたの回答
tips
プレビュー
