The twelve homework of computational physics of 5.1

孔德锋 2014301020081

---

摘要

因为我的spyer无法汇出三维图片（无函数包），因此此次作业采用matlab进行操作，汇出figure5.4的V的三维分布图

背景

使用书中介绍的the Jacobi method
对于二维空间使用：

$$V(i,j)=\frac{1}{4}(V(i+1,j)+V(i-1,j)+V(i,j+1)+V(i,j-1))$$
first：initial-V
second：update-V
third：laplace-calculate

正文及结论

利用matlab语言，可分成三步：
1. Laplace_prism
2. Initialise_prism
3. Update_prism

[V] =Initialise_prism;
% run update routine and estimate convergence
[V_new, delta_V_new]=Update_prism(V);
loops=0;
%Initialise loop counter
% While we have not met the convergence criterion and the number of loops is <10 so that we give the relaxation
% algorithm time to converge
while (delta_V_new>4e-5 | loops < 20);
    loops=loops+1;
    [V_new, delta_V_new]=Update_prism(V_new);
    mesh (V_new,'FaceColor','interp');
    title('Potential Surface');
    axis([0 20 0 20 0 1]);
    drawnow;
    pause(0.5);
end;

% This function creates the intial voltage array V
function[V] =Initialise_prism;
clear;
V = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
];

%This function creates the Update_prism array V
function [ V_new, delta_V_new ] = Update_prism(V);
row_size = size(V,1);
column_size=size(V,2);
V_new=V;
delta_V_new=0;
for j =2:column_size-1;
    for i=2:row_size -1;
% Do not update V in metal bar
        if V(i,j) ~=1;
% If the value of V is not =1, calculate V_new and
% delta_V_new to test for convergence
            V_new(i,j) = (V(i-1,j)+V(i+1,j)+V(i,j-1)+V(i,j+1))/4;
            delta_V_new=delta_V_new+abs(V_new(i,j)-V(i,j))
        else
% otherwise, leave value unchanged
            V_new(i,j)=V(i,j);
        end;
    end;
end;

感谢

Computational Physics using MATLAB®