Ta sẽ tính từ trường \(\vec{B}\) tại các điểm trong không gian giới hạn bởi \(\mathrm{xmin}\), \(\mathrm{xmax}\), \(\mathrm{ymin}\), \(\mathrm{ymax}\), \(\mathrm{zmin}\) và \(\mathrm{zmax}\):
1 2 3 4 5 6 7 8 9 |
Dây dẫn mang dòng điện \(I=1\mathrm{A}\) cần được mô tả sao cho sợi dây “trông có vẻ” dài vô hạn, mặc dù trên thực tế tính toán, ta chỉ có thể làm việc với đoạn dây hữu hạn. Để làm được điều đó, chỉ cần cho đoạn dây hữu hạn của chúng ta dài hơn vùng không gian được khảo sát một chút, vượt quá cả \(\mathrm{xmin}\) và \(\mathrm{xmax}\). Trong Matlab khai báo này thực thi qua lệnh:
1 2 3 4 |
Từ trường \(\vec{B}\) tại mỗi điểm trong không gian là sự chồng chập từ trường tạo ra bởi từng đoạn dây riêng rẽ, theo định luật Biot-Savart:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | Bx = zeros(size(X)); By = zeros(size(Y)); Bz = zeros(size(Z)); kB = mu0/4/pi*I; for i_wire = 1:N dx_wire = x_wire(i_wire+1)-x_wire(i_wire); dy_wire = y_wire(i_wire+1)-y_wire(i_wire); dz_wire = z_wire(i_wire+1)-z_wire(i_wire); Delta_x = X-x_wire(i_wire); Delta_y = Y-y_wire(i_wire); Delta_z = Z-z_wire(i_wire); r = sqrt(Delta_x.^2+Delta_y.^2+Delta_z.^2); Bx = Bx+kB*(dy_wire*Delta_z-dz_wire*Delta_y)./r.^3; % x-component of vector field By = By+kB*(dz_wire*Delta_x-dx_wire*Delta_z)./r.^3; % y-component of vector field Bz = Bz+kB*(dx_wire*Delta_y-dy_wire*Delta_x)./r.^3; % z-component of vector field end |

Code chương trình Matlab
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 | function M_field_straight_wire % Author: Tran Hai Cat % Lecturer in Physics, HCM University of Technology and Education % - Dai hoc Su pham Ky thuat Tp. Ho Chi Minh % Created: 2019.05.02 clc; clear variables; close all; %% CONSTANTS mu0 = 4*pi*1e-7; %% PARAMETERS I = 1; % Ampere x_wire = linspace(-2,4,50); y_wire = zeros(size(x_wire)); z_wire = zeros(size(x_wire)); streak_arrow = 1; % 1-MUI TEN CONG, 0-MUI TEN THANG Nx = 10; Ny = 6; Nz = 6; xmin = -1; xmax = 3; ymin = -1; ymax = 1; zmin = -1; zmax = 1; %% CALCULATION N = length(x_wire)-1; x = linspace(xmin,xmax,Nx); y = linspace(ymin,ymax,Ny); z = linspace(zmin,zmax,Nz); [X,Y,Z] = meshgrid(x,y,z); Bx = zeros(size(X)); By = zeros(size(Y)); Bz = zeros(size(Z)); kB = mu0/4/pi*I; for i_wire = 1:N dx_wire = x_wire(i_wire+1)-x_wire(i_wire); dy_wire = y_wire(i_wire+1)-y_wire(i_wire); dz_wire = z_wire(i_wire+1)-z_wire(i_wire); Delta_x = X-x_wire(i_wire); Delta_y = Y-y_wire(i_wire); Delta_z = Z-z_wire(i_wire); r = sqrt(Delta_x.^2+Delta_y.^2+Delta_z.^2); Bx = Bx+kB*(dy_wire*Delta_z-dz_wire*Delta_y)./r.^3; % x-component of vector field By = By+kB*(dz_wire*Delta_x-dx_wire*Delta_z)./r.^3; % y-component of vector field Bz = Bz+kB*(dx_wire*Delta_y-dy_wire*Delta_x)./r.^3; % z-component of vector field end %% FIGURES figure('name','Electric Field Rod','color','k','numbertitle','off'); if streak_arrow streakarrow3d(X,Y,Z,Bx,By,Bz,2,1); else quiver3(X,Y,Z,Bx,By,Bz,'AutoScaleFactor',3); end hold on % Draw the rod: plot3(x_wire,y_wire,z_wire,'linewidth',3); %% Axis properties axis equal axis([xmin xmax ymin ymax zmin zmax]); rotate3d on xlabel('X, m'); ylabel('Y, m'); zlabel('Z, m'); view(-20,30); box off grid off set(gca,'color','k','xcolor','w','ycolor','w','zcolor','w') function hh=streakarrow3d(X,Y,Z,U,V,W,np,AC) %% Bertrand Dano 03-09 % Copyright 1984-2009 The MathWorks, Inc. DX=abs(X(1,1,1)-X(1,2,1)); DY=abs(Y(1,1,1)-Y(2,1,1)); DZ=abs(Z(1,1,1)-Z(1,1,2)); DD=min([DX DY DZ]); ks=DD/100; % Size of the "dot" for the tuft graphs np=np*10; alpha = 3; % Size of arrow head relative to the length of the vector beta = .15; % Width of the base of the arrow head relative to the length XYZ=stream3(X,Y,Z,U,V,W,X,Y,Z); Vmag=sqrt(U.^2+V.^2+W.^2); Vmin=min(Vmag(:));Vmax=max(Vmag(:)); Vmag=Vmag(:); cmap=colormap; for k=1:length(XYZ) F=XYZ(k); [L M]=size(F{1}); if L<np F0{1}=F{1}(1:L,:); if L==1 F1{1}=F{1}(L,:); else F1{1}=F{1}(L-1:L,:); end else F0{1}=F{1}(1:np,:); F1{1}=F{1}(np-1:np,:); end P=F1{1}; if AC==1 vcol=floor((Vmag(k)-Vmin)./(Vmax-Vmin)*64); if vcol==0; vcol=1; end COL=[cmap(vcol,1) cmap(vcol,2) cmap(vcol,3)]; else COL='k'; end hh=streamline(F0); set(hh,'color',COL,'linewidth',.5); if L>1 x1=P(1,1); y1=P(1,2); z1=P(1,3); x2=P(2,1); y2=P(2,2); z2=P(2,3); u=x2-x1; v=y2-y1; w=z2-z1; uv=sqrt(u.*u + v.*v); xa1=x2+u-alpha*(u+beta*(v+eps)); xa2=x2+u-alpha*(u-beta*(v+eps)); xa3=x2+u-alpha*u; ya1=y2+v-alpha*(v-beta.*(u+eps)); ya2=y2+v-alpha*(v+beta.*(u+eps)); ya3=y2+v-alpha*v; za1=z2+w-alpha*w; za2=z2+w-alpha*(w+beta.*(uv+eps)); za3=z2+w-alpha*(w-beta.*(uv+eps)); plot3([x2 xa1 xa2 x2 xa3 xa3 x2],[y2 ya1 ya2 y2 ya3 ya3 y2],[z2 za1 za1 z2 za2 za3 z2],'color',COL); hold on end end axis image |