Scilabのエラー!二重振り子の数値解析で添字エラーが発生

このQ&Aのポイント
  • Scilabを使用して二重振り子の数値解析を行っていますが、エラーが発生しています。
  • オイラー法を使って運動方程式を解こうとしていますが、添字に誤りがあります。
  • ソースコードに問題がある可能性があります。20000行1列の行列を使用し、配列としてアクセスしています。どこに問題があるのか教えてください。
回答を見る
  • ベストアンサー

scilabのエラーに関して

よろしくお願いいたします。 Scilabを用いて二重振り子の数値解析を行なっています。 この先scilabが必要になるのでいろいろ使ってみているのですが、 うまくいきません。 オイラー法を用いて二重振り子の運動方程式を動かそうとしているのですが、 動きません。 ソースは //Matrix ccc1=zeros(20000,1); ccc2=zeros(20000,1); s1=zeros(20000,1); s2=zeros(20000,1); c1=zeros(20000,1); c2=zeros(20000,1); m1=0.10; //一重目質量 m2=0.08; //二重目質量 l1=0.4; //一本目長さ l2=0.3; //二本目長さ dt=0.05; //時間刻み幅 g=9.80665; //Gravitational acceleration c1(1)=%pi; c2(1)=%pi*0.5; s1(1)=0; //1. s2(1)=0; //2. //Eular method for t=0:1:100 z1=-m2*l1*s1(t+1)*s1(t+1)*sin(c1(t+1)-c2(t+1))*cos(c1(t+1)-c2(t+1)); z2=m2*g*cos(c1(t+1)-c2(t+1))*sin(c2(t+1)); z3=-m2*l2*s2(t+1)*s2(t+1)*sin(c1(t+1)-c2(t+1)); z4=(m1+m2)*g*sin(c1(t+1)); z5=l1*m1+l1*m2*sin(c1(t+1)-c2(t+1))*sin(c1(t+1)-c2(t+1)); ccc1=(z1+z2+z3+z4)/z5; ccc2(t+2)=ccc2(t+1)+dt*ccc2(t+1); s1(t+2)=s1(t+1)+dt*ccc1(t+1); c1(t+2)=c1(t+1)+dt*s1(t+1); z6=(sin(c1(t+1)-c2(t+1))*(m2*l2*s2(t+1)*s2(t+1)*cos(c1(t+1)-c2(t+1)))); z7=(sin(c1(t+1)-c2(t+1))*((m1+m2)*l1*s1(t+1)*s1(t+1))); z8=(m1+m2)*g*cos(c1(t+1)-c2(t+1))*sin(c1(t+1)); z9=(m1+m2)*g*sin(c2(t+1)); z10=(l2*m1+l2*sin(c1(t+1)-c2(t+1))*sin(c1(t+1)-c2(t+1))); ccc1=(z6+z7+z8+z9)/z10; ccc2(t+2)=ccc2(t+1)+dt*ccc2(t+1); s2(t+2)=s2(t+1)+dt*ccc2(t+1); c2(t+2)=c2(t+1)+dt*s2(t+1); end なのですが、 ccc2(t+2)=ccc2(t+1)+dt*ccc2(t+1); s1(t+2)=s1(t+1)+dt*ccc1(t+1); c1(t+2)=c1(t+1)+dt*s1(t+1); のあたりで、 「添字に誤りがあります」とエラーが出てしまいます。 20000行一列の行列を用意して、C言語のように配列を用いているのですが、、、。 ご教授よろしくお願いいたします。

質問者が選んだベストアンサー

  • ベストアンサー
  • kmee
  • ベストアンサー率55% (1857/3366)
回答No.1

ccc1=(z1+z2+z3+z4)/z5; としたものを ccc1(t+1); と使おうとしたからでしょう。

関連するQ&A

  • scilabのエラーについて

    よろしくお願いします。 scilabにてオーバーラップを掛けてFFT解析をしたく、いくつかのWebを参考にさせて頂きながら 試しているのですが、”添字に誤りがあります”とエラーになってしまいます。 ソースは clear all; stacksize('max') v0=read_csv('C:\Users\*****\Desktop\scilab\data7.csv',",") v1 = evstr(v0);//上記v0データは文字列認識の為数値データに変換 //窓関数 N1=256; win_l = window('hn',N1); win_l_minus = win_l * -1; //FFTを実行 N=size(v1,1);//データ長 N2=N1/4; //N/2で50%オーバーラップ L = floor(N/N2)-1;//全サンプルに対するFFTの必要回数  y = 0;//出力ベクトルyの初期値 past_tail=zeros(1,N2);//ハーフオーバーラップ加算信号の初期値 for k=1:L n =N2*(k-1);//FFTの開始地点を更新   v2=win_l*v1( n+1 : n+N1 );//観測信号に窓を掛ける v3=fft(v2,-1);//FFT y1=fft(v3,1);//IFFT y=(y,past_tail + y1(1:N2));       //今回のIFFTの前半を前回のIFFT後半に加算してから出力ベクトルyに追加 past_tail=y1(N2+1:N1);//今回のIFFT結果の後半を記憶 end fv1dBuV=zeros(N,2); fv1dBuV(1,1) = 0;//DC =0Hzを入力 fv1dBuV(1,2) = 20*log10(abs(y(1)/N)/10.^-6);//dBuVに変換 sr = 4000; //sr:サンプリングレート[sample/sec] dt = 1/sr; // dt: サンプリング間隔[sec] T = N*dt; // T: 測定時間[sec] df =1/T; //Hz df: 1/測定時間->最低周波数、周波数分解能[Hz] for i=2:N fv1dBuV(i,1)=(i-1)*df///10^6;//周波数をMHzで記録 fv1dBuV(i,2)=20*log10(abs((y(i))/(N/2))/10.^-6); //N/2の左半分を使用して算出 end になります。(この後はPlot処理) y=(y,past_tail + y1(1:N2))で上記エラーが指摘されます。 プログラミングは初心者なのでどのように修正したらよいか 教えていただけると助かります。 また、オーバーラップを掛けてのFFT解析で他の方法があれば 教えてください。

  • ラグランジュ方程式

    重ね重ねですが,ラグランジュ方程式に関して,本に以下のような式が例題として載っていました。微分するなかで,なぜそのような式に展開されるのか理解できません。よろしければ教えていただければと思います。 ∂L/∂v=(m1lG1^2+m2(l1^2+lG2^2)+2m2l1lG2cosθ2+I1+I2)v1+(lG2^2+m2l1lG2cosθ2+I2) という式をtで微分すると d(∂L/∂v)/dt=(m1lG1^2+m2(l1^2+lG2^2)+2m2l1lG2cosθ2+I1+I2)a1-2m2l1lG2v1v2sinθ2+(m2lG2^2+m2l1lG2cosθ2+I2)a2-m2l2lG2a2^2sinθ2 となると記載してありました。 この過程がよくわからないので,誰かご存知の方教えていただければと思います。

  • 複素解析の問題

    線分z=t*e^(π/4*i) (0≦t≦r)にそった ∫_C e^(-z^2)dzの積分の実部を、cos(t^2)とsin(t^2)を使って表せ この問題の答えは 1/√2(∫[0,r] cos(t^2)+sin(t^2) dt) で合っていますか?

  • 解析力学の2重振り子の問題について

    添付写真のような2重振り子の問題について考えているのですが ポテンシャルUはθ1=θ2=0のときを0とすると U=m1gl1(1-cosθ1)+m2g{l1(1-cosθ1)+l2(1-cosθ2)} となるそうです。 ここで困っているのですが、なぜ1-cosθiとなるのでしょうか? 私の考えでは U=m1gl1cosθ1+m2g(l1cosθ1+l2cosθ2)になると思います。 物理屋さんどなたかご教授おねがいします。

  • 定積分の結果に辿り着かせて下さい

    定積分の結果に辿り着かせて下さい。 概要 f_p(t) = { sin^2 t (0<=t<π) = { 0 (π<=t<2π) のフーリエ級数の f_s(t) = c_0 + Σ[n=1,∞] a_n cos nt + Σ[n=1,∞] b_n sin nt のb_nの定数を求めようとしています。 (本からの抜粋) b_m = (1/π) ∫[0,π] (sin t)^2 * sin mt dt (m = 1,2,3,...) の被積分関数は (sin t)^2 * sin mt = (1/2) sin mt - (1/4) {sin (2+m)t - sin (2-m)t} と書き直せるので、定積分することにより { b_m = 0 (mが偶数のとき) { b_m = - 4 / (πm(m^2 - 4)) (mが奇数のとき) が得られる。 …とありますが、計算間違いをしているのか、後一歩で辿り着けません。 私の計算 (1/2) ∫[0,π] sin mt dt - (1/4) {∫[0,π] sin (2+m)t dt - ∫[0,π] sin (2-m)t dt} =(1/2) [(1/m) (-cos mt)][0,π] - (1/4) { [(1/(2+m)) (-cos (2+m)t)][0,π] - (1/(2-m)) (-cos (2-m)t)][0,π] } =(-1/2m) [cos mt][0,π] - (1/4) { (-1/(2+m)) [cos (2+m)t][0,π] + (1/(2-m)) [cos (2-m)t][0,π] } =(-1/2m) [cos πm - cos 0] - (1/4) { (-1/(2+m)) [cos (2π+πm) - cos 0] + (1/(2-m)) [cos (2π-πm) - cos 0] } =(-1/2m) [cos πm - cos 0] - (1/4) { (-1/(2+m)) [cos πm - cos 0] + (1/(2-m)) [cos πm - cos 0] } =(-1/2m) [(-1)^(2m-1) - 1] - (1/4) { (-1/(2+m)) [(-1)^(2m-1) - 1] + (1/(2-m)) [(-1)^(2m-1) - 1] } =[ (-1/2m) - (1/4) { (-1/(2+m)) + (1/(2-m)) } ] * [(-1)^(2m-1) - 1] =[ (-1/2m) - (1/4) { (1/(2-m)) - (1/(2+m)) } ] * [(-1)^(2m-1) - 1] =[ (-1/2m) - (1/4) { (2+m)/{(2-m)(2+m)} - (2-m)/{(2-m)(2+m)} } ] * [(-1)^(2m-1) - 1] =[ (-1/2m) - (1/4) { (2+m-2+m)/(2-m)(2+m) } ] * [(-1)^(2m-1) - 1] =[ (-1/2m) - (1/4) { (2m)/(4-m^2) } ] * [(-1)^(2m-1) - 1] …とりあえず、ここまで合っていますでしょうか? 合っていたら、 { b_m = 0 (mが偶数のとき) { b_m = - 4 / (πm(m^2 - 4)) (mが奇数のとき) まで導いて下さい。 途中で計算間違いがあったらご指摘下さい。 お願いします。

  • 絶対値付き三角関数の積分、ラプラス変換の問題

    積分∫| cos(t) |e^-st dt を求めよ.という問題で 自分なりに計算してみたところ、 ∫| cos(t) |e^-st dt   (範囲は0~π) = ∫cos(t)e^-st dt - ∫cos(t)e^-st dt  (範囲は0~π/2、π/2~π) = [e^-st × (-scos(t) + sin(t) ) / s^2 + 1 ] - [e^-st × (-scos(t) + sin(t) ) / s^2 + 1] (範囲は0~π/2、π/2~π) = (2e^-πs/2 / s^2 + 1 ) - ( s / s^2 + 1 ) - (se^-πs / s^2 + 1) となりました。 その次の問題で、 |cos(t)| = |cos(t+π)|を用いて、ラプラス変換 L[ |cos(t)| ] = ∫|cos(t)|e^-st dt (範囲は0~∞)を計算せよ という問題があり、こちらの方は手も足も出ない状態です。。。 まず、前半の計算の仕方が合っているのかとその次の問題の解法をお伺いしたいです。 大変見にくいとは思いますが、どうかよろしくお願いします。

  • 積分の計算

    ∫1/√(x^2+1)dxをもとめよ。 x=tanθとおくと、dx=dθ/cos^2θ 与式=∫(dθ/cosθ)=∫cosθ/(1-sin^2θ)dθ sinθ=tとおくと、cosθdθ=dtより、 与式=∫dt/(1-t^2) =1/2((1/1-t)+(1/1+t))dt =1/2(-logI1-tI+logI1+tI)+C(絶対値) =1/2log{(1+t)/(1-t)}+C =1/2log{(1+sinθ)/(1-sinθ)}+C =1/2log{(1+sinθ)^2/cos^2θ}+C =log(1+sinθ/cosθ)+C とやって、tanθ=xを使って復元できなくなりました。 助けてください

  • 微分方程式

    L・d^2Q/dt^2+R・dq/dt+Q/C=V0sin(ωt) この解をQ(t)=Q0sin(ωt+θ)→I(t)=dq(t)/dt=Q0ωcos(ωt+θ) とおいて解くとtanθ=R/ω(L-C)となるのですがこれはどういう解法でこうなったのかがわかりません>< 分かる方は解法を教えて下さい。

  • RLC回路について

    RLC回路について 交流起電力 V=V?cosωtで L × dI/dt + RI +1/C × Q = V?cosωt 両辺をtで微分して dq/dt = I より L × d2I/dt2 + R × dIdt + 1/C×I = -ωV?sinωt この方程式でsinωtをどうにかオイラーの公式を使い実数で表して 特殊解 I = I?cos(ωt - δ) ただし I? = V?/{R? + (Lω - 1/Cω)?}?/?     tanδ = (Lω - 1/Cω)/R を求めるまでの過程が調べたのですがわかりません。 どなたかアドバイスお願いします。

  • 曲線の長さを求める問題が出来ません

    x = -1/3 (cos t)^3 y = sin t - 1/3 (sin t)^3 (0 ≦ t ≦ π/3) l = ∫[0, π/3]√(1- 2(cos t)^2 + 2(cos t)^4) dt 積分が出来ません。 どうやりますか?

専門家に質問してみよう