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
|
N=100;
Deltatheta=pi/5;
Deltap=1;
theta=linspace(-2*pi,2*pi,N);
imp=linspace(-3,3,N);
[Theta Imp]=meshgrid(theta, imp);
f=zeros(N);
%//////////////////////////////////////////
%conditions initiales/////////////
f((3*N)/10:(7*N)/10,N/6:(5*N)/6)=60/(N*N);
Theta0=Theta;
Imp0=Imp;
%/////////////////////////////////
%surface(f)
deltat=0.01;
%iteration sur le temps/////////////
for i=1:100
M=trapz(imp,trapz(theta,cos(Theta0).*f)); %calcul des intégrales
M2=trapz(imp,trapz(theta,sin(Theta0).*f));
DV=sin(Theta0).*M - cos(Theta0).*M2; % calcul de la dérivée du potentiel
Theta0= Theta0 + deltat.*Imp0;
Imp0= Imp0 - deltat.*DV;
f1=griddata(Theta0(:),Imp0(:),f(:),Theta,Imp,'linear');
f=f1;
end
%//////////////////////////////////////
surface(f) |
Partager