四 计算平面上N个电荷之间的库伦引力
1 建模：
由库仑定律： F?q1q2/4??0r2 Fx?q1q2(x2?x1)/4??0r3 3F?qq(y?y)/4??ry12210
r?(x2?x1)2?(y2?y1)2
先输入电荷的数目，各电荷的坐标及电荷量，再选一个电荷，求其它电荷对它的作用力，叠加求合力。再选下一个电荷，依次类推。
Matlab程序：
clear all;
N = input('ê?è?μ?oéêy??N=:'); for ic = 1:N fprintf('----/n??μ?oé#%g\\n',ic); rc = input('ê?è?μ?oé????[x,y]￡¨?×￡?:'); x(ic) = rc(1); y(ic) = rc(2);
q(ic) = input('ê?è?μ?oéá?￡¨?a??￡?￡o'); end
E0 = 8.85e-12; C0 = 1/(4*pi*E0); for ic = 1:N Fx = 0.0;Fy = 0.0; for jc = 1:N
if(ic ~= jc)
xij = x(ic)-x(jc);yij = y(ic)-y(jc); Rij = sqrt(xij^2+yij^2); Fx=Fx+C0*q(ic)*q(jc)*xij/Rij^3; Fy=Fy+C0*q(ic)*q(jc)*yij/Rij^3; end end
fprintf('???üμ?oé×÷ó??úμ?oé#%gé?μ?o?á|?a￡o\\n',ic); fprintf('x-·?á?:%e\\n',Fx); fprintf('y-·?á?:%e\\n',Fy); end
五 有限差分法处理电磁场问题
Matlab程序:
m=40 for k=1:m for j=1:m if k==1 V(j,k)=1;
elseif((j==1)|(j==m)|(k==m)) V(j,k)=0; else
V(j,k)=0.5; end end end cha=0.01; delta=0; n=0;
while(1) n=n+1; for k=2:m-1 for j=2:m-1
Vnew(j,k)=1/4*(V(j+1,k)+V(j-1,k)+V(j,k+1)+V(j,k-1)); d=abs((Vnew(j,k)-V(j,k))/V(j,k)); if d>delta delta=d; end
V(j,k)=Vnew(j,k); end end
if delta>cha break; end if(n>100) break; end delta=0.; end
结果为：