前提・実現したいこと
c言語でルンゲクッタ法を用いて出したデータをプロットしたいです。
プログラム内でファイルに書き出すやり方です。
発生している問題・エラーメッセージ
gnuplotでプロットしようとしたら以下のエラーが発生しました。
warning: Skipping data file with no valid points ^ x range is invalid
該当のソースコード
#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"); if(fp==NULL){ printf("cannot open\n"); fclose(fp); } } 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); }
データを貼り付けてプロットする方法も教えていただきたいです。
補足情報(FW/ツールのバージョンなど)
言語:C Xcodeを使っています。
プロットはgnuplotを使って行いたいです。
10個のブロックをバネで繋げたモデルで、走る際の荷重の掛け方による摩擦の変化を見るプログラムです。
つま先、足裏全体、踵の3パターンをθなどの値を変えて調べたいと思っています。
お願い致します。
回答2件
あなたの回答
tips
プレビュー