一次元静電粒子コードについて・・ -これは、VBで作ってみたPICコードな-

一次元静電粒子コードについて・・

一度物理学のところで質問したのですが、誰も答えてくれなかったので(ずうずうしいとは思いますが)今度はこちらで質問させていただきます。これは、VBで作ってみたPICコードというものなのですが(電界と粒子の運動の関係のコードです)、このコードでは右側の図にSin(正弦)波形が出力されます。この波形を変えたいのですが、どこを変えればいいのでしょう？私は、ポテンシャルの初期値を変えればいいと思い何度もチャレンジしましたが、うまくいきません。　どなたか分かる方教えてください!! （つづきがあります。） '一次元静電粒子コード M = 3:IM = 2 ^ (M + 1) + 1:IM2 = IM * 2 - 1 Dim RO(18):Dim PHI(18):Dim A(18, 4):Dim E(18):Dim X(161):Dim V(161) ' X:電子の位置（Ｘ座標） V:電子の速度 RO:電荷密度 ' A(I,J):ポアソン・マトリックスエレメント PHI:静電ポテンシャル ' E:静電界（X方向） NP:1セル当りの粒子数 XMAX:全メッシュ数 ' KK:最大計算回数 DT:タイムステップ XMAX:X方向の最大値点 ' DX: メッシュ間隔 NP = 10:NPT = NP * (IM - 1):PW = 1 / NP XMAX = IM - 1:DX = XMAX / (IM - 1):DXP = DX / NP:DT = 0.1 ' 粒子の初期設定 For IP = 1 To NPT X(IP) = DXP * (IP - 0.5):V(IP) = 0 Next IP ' ポテンシャルの初期値 For I = 1 To IM PHI(I) = 5! * Sin(2 * 3.14159 * (I - 1) / (IM - 1)) Next I ' 計算開始 For KK = 1 To 120 ' 電場の計算 For I = 1 To IM - 1 E(I) = -(PHI(I + 1) - PHI(I)) / DX Next I E(IM) = E(1):E(0) = E(IM - 1):E(IM + 1) = E(2) つづく・・・

一次元静電粒子コードについて・・(つづき１)

' 粒子の運動 For IP = 1 To NPT XNM = X(IP) / DX + 0.5:IXP = Int(XNM):IXM1 = IXP - 1 If IXM1 <= 0 Then IXM1 = IM - 1 - IXM1: B = XNM - IXP If IXP = 0 Then IXP = IM - 1:EP = E(IXP) * B + E(IXM1) * (1 - B) V(IP) = V(IP)-EP * DT: X(IP) = X(IP) + V(IP) * DT If X(IP) > XMAX Then X(IP) = X(IP) - XMAX If X(IP) < 0 Then X(IP) = XMAX + X(IP) Next IP ' 電荷密度の計算 For I = 1 To IM - 1 RO(I) = -1 Next I For IP = 1 To NPT XNM = X(IP) / DX:IXP = Int(XNM) B = XNM - IXP :If IXP = 0 Then IXP = IM - 1 IXP1 = IXP + 1 :If IXP1 = IM Then IXP1 = 1 RO(IXP) = RO(IXP) + PW * (1 - B):RO(IXP1) = RO(IXP1) + PW * B Next IP RO(IM) = RO(1) For I = 1 To IM - 1 RO(I) = RO(I) * DX * DX Next I ' 静電ポテンシャルの計算 For I = 2 To IM - 1 A(I, 1) = 1:A(I, 2) = -2:A(I, 3) = 1:A(I, 4) = RO(I) Next I A(1, 4) = 0 :A(1, 2) = 1:A(IM, 4) = 0:A(IM, 2) = 1:BO = -2 ' -------ガウスの消去法------- '前進消去法 NN = 1:N = IM:NNP1 = NN + 1:NN2 = 2 * NNP1 For I = 1 To N - 1 For J = NN2 To NNP1 Step -1 A(I, J) = A(I, J) / A(I, NNP1) Next J IPNN = I + NN:If IPNN > N Then IPNN = N For L = I + 1 To IPNN For K = NNP1 + I - L + 1 To 2 * NN + 1 + I - L A(L, K) = A(L, K) - A(I, K + L - I) * A(L, NNP1 + I - L) Next K A(L, NN2) = A(L, NN2) - A(I, NN2) * A(L, NNP1 + I - L) Next L Next I つづく関連URL:http://www.okweb.ne.jp/kotaeru.php3?q=449543

一次元静電粒子コードについて・・(つづき１)

' 電荷密度の計算 For I = 1 To IM - 1 RO(I) = -1 Next I For IP = 1 To NPT XNM = X(IP) / DX:IXP = Int(XNM) B = XNM - IXP :If IXP = 0 Then IXP = IM - 1 IXP1 = IXP + 1 :If IXP1 = IM Then IXP1 = 1 RO(IXP) = RO(IXP) + PW * (1 - B):RO(IXP1) = RO(IXP1) + PW * B Next IP RO(IM) = RO(1) For I = 1 To IM - 1 RO(I) = RO(I) * DX * DX Next I ' 静電ポテンシャルの計算 For I = 2 To IM - 1 A(I, 1) = 1:A(I, 2) = -2:A(I, 3) = 1:A(I, 4) = RO(I) Next I A(1, 4) = 0 :A(1, 2) = 1:A(IM, 4) = 0:A(IM, 2) = 1:BO = -2 ' -------ガウスの消去法------- NN = 1:N = IM:NNP1 = NN + 1:NN2 = 2 * NNP1 For I = 1 To N - 1 For J = NN2 To NNP1 Step -1 A(I, J) = A(I, J) / A(I, NNP1) Next J IPNN = I + NN:If IPNN > N Then IPNN = N For L = I + 1 To IPNN For K = NNP1 + I - L + 1 To 2 * NN + 1 + I - L A(L, K) = A(L, K) - A(I, K + L - I) * A(L, NNP1 + I - L) Next K A(L, NN2) = A(L, NN2) - A(I, NN2) * A(L, NNP1 + I - L) Next L Next A(N, NN2) = A(N, NN2) / A(N, NNP1) For I = N To 2 Step -1 IMNN = I - NN:If IMNN < 1 Then IMNN = 1 For L = I - 1 To IMNN Step -1 A(L, NN2) = A(L, NN2) - A(I, NN2) * A(L, NNP1 + I - L) Next L Next I A(1, NN2) = A(1, NN2) / A(1, NNP1) For I = 1 To IM PHI(I) = A(I, 4) Next I つづく関連URL:http://www.okweb.ne.jp/kotaeru.php3?q=448303

一次元静電粒子コードについて・・(つづき２)

'後退代入法 A(N, NN2) = A(N, NN2) / A(N, NNP1) For I = N To 2 Step -1 IMNN = I - NN:If IMNN < 1 Then IMNN = 1 For L = I - 1 To IMNN Step -1 A(L, NN2) = A(L, NN2) - A(I, NN2) * A(L, NNP1 + I - L) Next L Next I A(1,NN2)=A(1,NN2)/A(1,NNP1) 'グラフ Picture1.BackColor = RGB(255, 255, 255) Picture1.Line (50, 100)-(50, 300):Picture1.Line (50, 200)-(250, 200) Picture1.Line (350, 100)-(350, 300) Picture1.Line (350, 200)-(550, 200) '電子の位置AXに対する速度AYの変化　　(LEFT SIDE) For I = 1 To NPT 'NPT=NP*(IM-1):全メッシュ数=160個 AX = X(I) * 200 / (IM - 1) + 50 :AY = 200 - V(I) * 100 Picture1.Circle (AX, AY), 1, vbBlue Next I '電界の初期値 AXO = 350:AYO = 200 + E(1) * 100 Picture1.Circle (AXO, AYO), 10, vbGreen '電子の位置（格子点）に対する電界の値 (RIGHT SIDE) For I = 2 To IM - 1 AX = 350 + (I - 1) * 200 / (IM - 1):AY = 200 + E(I) * 100 Picture1.Circle (AX, AY), 1, vbGreen Picture1.Line (AXO, AYO)-(AX, AY), vbRed;AXO = AX :AYO = AY Next I Next KK End Sub 以上です。長ったらしく、また見にくくて申し訳ありません。　分かる方がいましたら、ぜひ教えてください。関連URL:http://www.okweb.ne.jp/kotaeru.php3?q=449543 関連URL:http://www.okweb.ne.jp/kotaeru.php3?q=449550

C言語　初心者です

Cで計算して出た結果をエクセルに出力するやり方を教えて下さい。 #include<stdio.h> #include<math.h> double V(double x) { double pot=0.0; if (x >-1 && x<1) { pot=-2.0; } return(pot); } double F(double psi, double E,double x) { return(2.*(V(x)-E)*psi); } int main() { double E; int i,j; double psi,psi1,phi,phi1,x,dx; psi=0.01; phi=0.01; scanf("%1f",& E); x=-5.; i=5; dx=pow(10,-i)*10.0; for(j=0;j<pow(10,i);j++) { psi1=psi+phi*dx; phi1=phi+F(psi,E,x)*dx; x=x+dx; psi=psi1; phi=phi1; if((j+1)/5000*5000==j+1) { printf("x=%f phi= %f psi= %f\n",x,phi,psi); } } }

ベーシックをVBに訳してください。

print "time=",kk　　　　　　　 for i=1 to im print using 　"i=## 　ro=###.###~~~~　　phi=###.###~~~~　E=###.###~~~~" ;i;ro(i);phi(i);e(i) next i どなたかこのベーシックで書かれたコードを、VBに直していただけないでしょうか？　お願いします。

2次元極座標の速度

2次元平面内を運動する物体がある。この物体の運動をデカルト座標を用いて表すと、その速度は v(t)={dx(t)/dt}i+{dy(t)/dt}j　　(v(t)は速度ベクトル) と表せる。但し、i,jはx軸y軸それぞれの単位ベクトルを表すものとする。同様のことを２次元極座標を使って表すと、動径方向の単位ベクトルをer(t)、角度方向の単位ベクトルをeθ(t)として、 v(t)={dr(t)/dt}er(t)+r{dθ(t)/dt}eθ(t) と書けることを示せ。但し、r(t)=√{x(t)^2+y(t)^2}とし、角度θ(t)は x軸から測った角度とする。(x(t)=r(t)cosθ(t),y(t)=r(t)sinθ(t)) という問題ですが、 (解)(長くなるので途中の式は省きます。) dx(t)/dt={dr(t)/dt}cosθ(t)-r(t){dθ(t)/dt}sinθ(t) dy(t)/dt={dr(t)/dt}sinθ(t)+r(t){dθ(t)/dt}cosθ(t) とそれぞれ求め、v(t)={dx(t)/dt}i+{dy(t)/dt}jに代入すると v(t)={dr(t)/dt}{cosθ(t)i+sinθ(t)j}+r(t){dθ(t)/dt}{-sinθ(t)i+cosθ(t)j}・・・アと表せる。 er(t)=Ai+Bj eθ(t)=Ci+DjとおいてA,B,C,Dをもとめると er(t)=cosθ(t)i+sinθ(t)j・・・イ eθ(t)=-sinθ(t)i+cosθ(t)j・・・ウアイウよりv(t)={dr(t)/dt}er(t)+r{dθ(t)/dt}eθ(t)と表せる。という解き方をしたんですが、適切ですか？お願いします。もっといい解き方があれば教えてください。

水中の微粒子分布はポアソン分布になるのでしょうか？

「100mlの水に500個の微粒子を入れ、均一になるように良く撹拌してあります。ここから10ｍｌすくい取ったとき、ｘ個の微粒子が存在する確率を求めたい」という場合、ポアソン分布になっているのでしょうか。 100ml中に500個では微粒子数が多すぎてポアソン分布になっていないような気がするのですが、10^5μl中に500個あると考えるとポアソン分布でいいような気もします。私は、ポアソン分布の確率関数 f(x)=e^(-λ)*λ^x/x! において、 n:サンプル量(μl)　p：微粒子濃度(個/μl)　λ=np とし、n=10^4、 p=5*10^(-3)、 λ=50 より f(x)=e^(-50)*50^x/x! と考えたのですが、合っているでしょうか？容量の単位を変えると微粒子濃度が大きくなったり小さくなったり感じられ、ポアソン分布の適用基準がわかりません。本などで調べたのですが類似の例がなく、良くわかりません。宜しくお願いいたします。

Scilabを使ったジュリア集合の描画プログラム

今、Scilabを使用してジュリア集合（充填およびそれ以外を含む集合）を描画するプログラムを書いています。以前書いたC言語のプログラムをもとに書いているのですが、正確に描画できません。どうしたらよいでしょうか。教えてください。与える条件は右上と左下の座標（複素数形式）と定数Cの値です。以下に掲載したのが製作したプログラムです。よろしくお願い致します。 //描画エリアの右上と左下の座標を複素数で設定する Z =[ -1.2-1.2*%i;1.2+1.2*%i]; //複素定数(C)を設定する C=0+0*%i; //描画エリアのx座標とy座標の各最小値と最大値を計算する。 xmin = min(real(Z)); xmax = max(real(Z)); ymin = min(imag(Z)); ymax = max(imag(Z)); Cr = real(C); Ci = imag(C); //描画点数を800×800に設定する。 N = 800; //各増分を計算する。 dx = (xmax-xmin)/(N-1); dy = (ymax-ymin)/(1-N); //プロットデータを"0"で初期化 map=zeros(N,N); //ジュリア集合の描画 i=1; for X=xmin:dx:xmax j=1; for Y=ymax:dy:ymin for k=1:30 x = X ^ 2 - Y ^ 2 + Cr; y = 2 * X * Y +Ci; if x^2 +y^2 > 4 then break; end map(j,i)=k; X=x; Y=y; end j=j+1; end i=i+1; end //プロットするための設定 Re = xmin:dx:xmax; Im = ymax:dy:ymin; clf(0); Sgrayplot(Re,Im,map');

合成関数の積分方法

久しぶりの積分でかなり忘れてるのでよろしくお願いします． I = ∮e^(-x^2)dx,D : 0≦x≦1についてですが， -1/2［e^(x^2)］D=(1-e)/2ってできますっけ？最初置換積分で解こうとして， x^2=tと置き，x=√t(∵x≧0), dx = 1／(2√t)dt, 0≦t≦1より I = ∮(e^(-t)・1/2√t)dtとなったんですが，部分積分法で解けませんでした．

置換積分の途中計算がわかりません

ー教科書－－－－－－－－－－－－－－－－－－－ (e^x)+1=tとおいて置換積分すると ∫((e^x)+1)' log((e^x)+1)dx =∫(logt)(dt/dx)dx ー－－－－－－－－－－－－－－－－－－－とありましたが、 (e^x)+1=tを全微分すると (e^x)dx=dtより dx=(1/(t-1))dtとなるためー－－－－－－－－－－－－－－－－－－－－－ ∫((e^x)+1)' log((e^x)+1)dx =∫t' (logt)(1/(t-1))dt =∫(logt)(1/(t-1))dt ー－－－－－－－－－－－－－－－－－－－－－ではないのですか？どこかで私の計算が間違っているのだと思われます。よろしくお願いします。

微分方程式の一般解

(i) d^2x/dt^2+3dx/dt+2x=t (ii)d^2x/dt^2-2dx/dt=6 途中式なども詳しく教えていただけるとうれしいです。ついでに私は最初にx=e^λtと仮定して、代入して解こうと思ったのですがこのやり方は間違ってますか？どの微分方程式にもこの方法はつかえるのですか？よろしくおねがいします。

偏微分方程式、差分法

Fitzhugh-Nagumo方程式 dU/dt = d/dx(dU/dx) + (a-U)(U-1)U-V (d/dx(dU/dx)はUのxに対する二回微分) dV/dt =e(bU-gV) a=0.1, b=0.5, g=1.0, e=0.01 初期条件 (U,V)=(1,0) if x=0 (U,V)=(0,0) if x>0 境界条件 dU/dx=0 at x=0 dV/dx=0 at x=0 を差分して陽解法で解くと添付の答えになるらしいのですが、自分で解いたところ添付の結果のようになりました。答えと一致しないため、プログラム上で何を間違えているのかご指摘頂けると助かります。 #include <stdio.h> #include <stdlib.h> //exsit()で必要 #include <math.h> int main(){ double a,b,g,e; double dt,dx; int x,t; double **u,**v; int i,j,k; a=0.1;b=0.5;g=1.0;e=0.01; x=100;t=5000.; dt=0.001;dx=0.1; FILE *fp; if((fp = fopen("ResultNagumo.txt","w"))==NULL){ printf("Can't open file\n"); exit(2); } u = (double**)malloc(sizeof(double*)*t); v = (double**)malloc(sizeof(double*)*t); for(i=0;i<t;i++){ u[i]=(double*)malloc(sizeof(double)*x); v[i]=(double*)malloc(sizeof(double)*x); } //初期値 u[0][0]=1.0; v[0][0]=0.0; for(i=1;i<x;i++){ u[0][i]=0.0; v[0][i]=0.0; } 　　 //差分計算 for(i=0;i<t-1;i++){ for(j=1;j<x-1;j++){ u[i+1][j] = u[i][j] + dt*((u[i][j+1]-2*u[i][j]+u[i][j-1])/pow(dx,2)+(a-u[i][j])*(u[i][j]-1)*u[i][j]-v[i][j]); v[i+1][j] = v[i][j] + dt*(e*(b*u[i][j]-g*v[i][j])); } //境界条件 u[i+1][0]=u[i+1][1]; v[i+1][0]=v[i+1][1]; u[i+1][x-1]=u[i+1][x-2]; v[i+1][x-1]=v[i+1][x-2]; } 　　//結果の出力 for(i=0;i<t;i++){ if((i%100)==0){ fprintf(fp,"%d\n",i); for(j=0;j<x;j++){ fprintf(fp,"%2.4e,",u[i][j]); } } } for(i=0;i<t;i++){ free(u[i]); free(v[i]); } free(u);free(v); fclose(fp); getchar(); return 0; }

微分の質問です。

途中式を書いてもらえると嬉しいです。よろしくお願いします線形方程式の初期値問題をラプラス変換してください。１ d^2x/dt^2-3(dx/dt)-10x=0, x(0)=0 dx/dt(0)=7 2 d^2x/dt^2-4(dx/dt)+4x=0, x(0)=1 dx/dt(0)=0 3 d^2x/dt^2-3(dx/dt)+2x=e^3t, x(0)=1 dx/dt(0)=0

微分の問題です

途中式を書いてもらえると嬉しいです。よろしくお願いします線形方程式の初期値問題をラプラス変換でといてください。１ d^2x/dt^2-3(dx/dt)-10x=0, x(0)=0 dx/dt(0)=7 2 d^2x/dt^2-4(dx/dt)+4x=0, x(0)=1 dx/dt(0)=0 3 d^2x/dt^2-3(dx/dt)+2x=e^3t, x(0)=1 dx/dt(0)=0

微分の問題です。

途中式を書いてもらえると嬉しいです。今日の試験で似たような問題が出るので急ぎでよろしくお願いします線形方程式の初期値問題をラプラス変換してください。１ d^2x/dt^2-3(dx/dt)-10x=0, x(0)=0 dx/dt(0)=7 2 d^2x/dt^2-4(dx/dt)+4x=0, x(0)=1 dx/dt(0)=0 3 d^2x/dt^2-3(dx/dt)+2x=e^3t, x(0)=1 dx/dt(0)=0

【問題】∫[0 to 1](e^x-1)/(e^(2x)+e^(-x)

【問題】∫[0 to 1](e^x-1)/(e^(2x)+e^(-x))dxを計算せよ。 e^x=tとおいて、e^x*dx=dt これより、∫[1 to e](t-1)/(t^3+1)dtとおけたのですが…こっからどうにもできません^^; どなたかよろしくお願いします!!!!

原始関数の問題の解き方

以下のように解いたのですが、解答に自信がありません。途中の式など、間違っていればご指摘のほどよろしくお願いします。次の原始関数を求めよ。 (1) ∫(x+1)^5 dx x+1=tとおく。 (dt/dx)=1より、dx=dt よって、∫(x+1)^5 dx=∫t^5 dt =(1/6)t^6+C =(1/6)(x+1)^6+C (Cは積分定数) (2) ∫e^(5x) dx 5x=tとおく (dt/dx)=5より、dx=(dt/5) =∫e^(t)(dt/5)+C =(1/5)e^(5x)+C (Cは積分定数) (3) ∫x/(x^2+1)^2 dx =∫{(x+1)-1}/(x^2+1)^2 dx =(1/2)∫{(2x+2)-2}/(x^2+1)^2 dx =(1/2)∫(x^2+1)'/(x^2+1)^2 dx =(1/2)log|(x^2+1)^2|+C (Cは積分定数) (4) ∫1/√(23-x^2) dx 公式 ∫1/√(a^2-x^2) dx=sin^(-1) x/√a+C (a>0)より =sin^(-1) x/√23 +C (Cは積分定数) ご指導、よろしくお願いします。

∫dx/√(e^x-1)の計算

∫dx/√(e^x-1)の計算なのですが、私はこう解きました。 t=√(e^x-1)とおいてdt=1/2tdx で、 ∫dx/√(e^x-1) =∫1/t・2tdt =∫2dt =2t =2√(e^x-1) としたのですが、答えは2tan^-1√(e^x-1)でした。どこが違うのか分かる方お答えお願いします。

積分　証明　問題

積分　証明　問題 ∫[0～π]（x・sinx）dxを求めよ。 I=∫[0～π]（x・sinx）dxとおく。 x=π-tとおくと、dx/dt=-1、積分範囲はπ～0 I=∫[π～0]（π-t）・sin（π-t）（-dt） =∫[0～π]（π-t）・sin（π-t）dt =∫[0～π]（π-t）・(sint)dt 2I=∫[0～π]（x・sinx）dx+∫[0～π]（π-x）・(sinx)dt 　=∫[0～π]πsinxdx 　=2π I=π 一点分からない点があります。 ∫[0～π]（π-t）・(sint)dt=∫[0～π]（π-x）・(sinx)dt について。単純にtをxに置き換えただけだと思いますが、 x=π-tと置換しているのに、t=xと同じ変数を使って再度置換して良いのでしょうか？以上、ご回答よろしくお願い致します。

一次元静電粒子コードについて・・

みんなの回答

関連するQ&A

一次元静電粒子コードについて・・

一次元静電粒子コードについて・・(つづき１)

一次元静電粒子コードについて・・(つづき１)

一次元静電粒子コードについて・・(つづき２)

C言語　初心者です

ベーシックをVBに訳してください。

2次元極座標の速度

水中の微粒子分布はポアソン分布になるのでしょうか？

Scilabを使ったジュリア集合の描画プログラム

合成関数の積分方法

置換積分の途中計算がわかりません

微分方程式の一般解

偏微分方程式、差分法

微分の質問です。

微分の問題です

微分の問題です。

【問題】∫[0 to 1](e^x-1)/(e^(2x)+e^(-x)

原始関数の問題の解き方

∫dx/√(e^x-1)の計算

積分　証明　問題

注目のQ&A

カテゴリ
一覧

専門家に質問してみよう
専門家登録

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

一次元静電粒子コードについて・・

みんなの回答

関連するQ&A

一次元静電粒子コードについて・・

一次元静電粒子コードについて・・(つづき１)

一次元静電粒子コードについて・・(つづき１)

一次元静電粒子コードについて・・(つづき２)

C言語 初心者です

ベーシックをVBに訳してください。

2次元極座標の速度

水中の微粒子分布はポアソン分布になるのでしょうか？

Scilabを使ったジュリア集合の描画プログラム

合成関数の積分方法

置換積分の途中計算がわかりません

微分方程式の一般解

偏微分方程式、差分法

微分の質問です。

微分の問題です

微分の問題です。

【問題】∫[0 to 1](e^x-1)/(e^(2x)+e^(-x)

原始関数の問題の解き方

∫dx/√(e^x-1)の計算

積分 証明 問題

注目のQ&A

カテゴリ 一覧

専門家に質問してみよう 専門家登録

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

C言語　初心者です

積分　証明　問題

カテゴリ
一覧

専門家に質問してみよう
専門家登録