% PROGRAMME D'ECRITURE POUR LA RUPTURE % declaration des variables globales global nc REPAR X Vrup BI1 BI2; % GENERATION DES CLASSES % volume minimum Vmin=1; % definition du nombre de classes nc=25; % raison geometrique pour la descretisation des classes s=1.5; % vecteur des volumes de discretisation X=zeros(nc,1); for i=1:nc X(i)=Vmin*s^(i-1); end % CALCULS PRELIMINAIRES % exposent pour loi de conservation e1=0; e2=1; % AGGLOMERATION % test pour la repartition pour agglomeration % initialisation des tableaux m=3*nc; POINTE=ones(nc,1); REPAR=zeros(nc,m,3); for i=1:nc for j=i:nc xf=X(i)+X(j); for k=j:nc-1 if (xf >= X(k)) & (xf < X(k+1)) REPAR(k,POINTE(k),1)=i; REPAR(k,POINTE(k),2)=j; REPAR(k,POINTE(k),3)=(xf^e1*X(k+1)^e2-xf^e2*X(k+1)^e1)/(X(k)^e1*X(k+1)^e2-X(k)^e2*X(k+1)^e1); POINTE(k)=POINTE(k)+1; REPAR(k+1,POINTE(k+1),1)=i; REPAR(k+1,POINTE(k+1),2)=j; REPAR(k+1,POINTE(k+1),3)=(xf^e1*X(k)^e2-xf^e2*X(k)^e1)/(X(k+1)^e1*X(k)^e2-X(k+1)^e2*X(k)^e1); POINTE(k+1)=POINTE(k+1)+1; end end if xf==X(nc) REPAR(nc,POINTE(nc),1)=i; REPAR(nc,POINTE(nc),2)=j; REPAR(nc,POINTE(nc),3)=1; POINTE(nc)=POINTE(nc)+1; end end end % RUPTURE % Vecteur des frequences de rupture Vrup=zeros(nc,1); for i=1:nc Vrup(i)=frup(X(i)); end % Evaluation des integrales I et II pour tous les couples BI1=zeros(nc,nc); BI2=zeros(nc,nc); % definition des 2 variables symboliques syms v1 real; syms v; % fonction de repartition frepar=2/v1; % fonction a integrer f=v^e1*frepar; g=v^e2*frepar; % Evaluatioms des integrales I for i=1:nc-1 for k=i+1:nc %v1=sym(X(k)); v1=X(k); I1= eval(int(f,v,X(i),X(i+1))); I2= eval(int(g,v,X(i),X(i+1))); BI1(i,k)=(X(i+1)^e2*I1-X(i+1)^e1*I2)/(X(i+1)^e2*X(i)^e1-X(i+1)^e1*X(i)^e2); end end % Evaluatioms des integrales II for i=2:nc for k=i:nc %v1=sym(X(k)); v1=X(k); I1= eval(int(f,v,X(i-1),X(i))); I2= eval(int(g,v,X(i-1),X(i))); BI2(i,k)=(X(i-1)^e2*I1-X(i-1)^e1*I2)/(X(i-1)^e2*X(i)^e1-X(i-1)^e1*X(i)^e2); end end % declaration des variables symboliques syms v1 real; syms v; % INITIALISATION %Distribution initiale %di=v*exp(-v); % Distribution initiale des classes N=zeros(nc,1); N(1)=1; % Pour tracage des courbes en surface %for i=1:nc-1 % N(i)=eval(int(di,v,X(i),X(i+1))); % S(i)=N(i)*pi*X(i)*X(i); %end %N(nc)=eval(int(di,v,X(nc),s*X(nc))); %S(nc)=N(nc)*pi*X(nc)*X(nc); % RESOLUTION tf=10; % nombre entier multiple de 10 %en nombre [x,y] =ode23s('equat',[0:tf],N); %en surface %[x,y] =ode23s('equat',[0:tf],S); sx=max(size(x)); t1=floor(2*sx/10); t2=floor(4*sx/10); t3=floor(6*sx/10); t4=floor(8*sx/10); t5=floor(10*sx/10); semilogx(X,y(t1,:)/sum(y(t1,:)),'r+-',X,y(t2,:)/sum(y(t2,:)),'go-',X,y(t3,:)/sum(y(t3,:)),'b*-',X,y(t4,:)/sum(y(t4,:)),'cx-',X,y(t5,:)/sum(y(t5,:)),'yO-'); %plot(X,y(5,:),'r+',X,y(10,:),'go',X,y(15,:),'b*',X,y(20,:),'cx',X,y(25,:),'k*'); %axis([0, 10, -inf, inf]); hold on semilogx(X,N,'m*-'); %plot(X,N,'-m'); %plot(X,y(2,:),X,y(3,:),X,y(5,:),X,y(7,:),X,y(10,:)) hold off