Отправлено: 31.05.10 22:26. Заголовок: Стабилизатор в R3

По многочисленным просьбам трудящихся ниже приведён код двух вариантов счёта стабилизирующего функционала в пространстве R3.
Первый-исходный, написанный Андреем Михайловичем Волковым. Второй мой. Суть у них абсолютно одинаковая, но реализация существенно отличается.
aq32.m (Волков edition)

 
 function[st]=aq32(n1,n2,n3) 
 p=ones(1,n1);q=1:1:n2;l=p;ll=q; 
 for i=2:n1 
    p=[p l*i];q=[q ll]; 
 end 
 p=[p;q]';nn=n1*n2;r=ones(nn,1); 
 pp=[p r]; 
 for i=2:n3 
    pp=[pp;p r*i]; 
 end 
 np=size(pp); 
 t0x=q0(n1);t0y=q0(n2);t0z=q0(n3); 
 %t1x=q1(n1);t1y=q1(n2);t1z=q1(n3); 
 t2x=q2(n1);t2y=q2(n2);t2z=q2(n3); 
 for i=1:n1*n2*n3 
    for j=1:n1*n2*n3 
       nn=[pp(i,:) pp(j,:)]; 
       d1=t2x(nn(1),nn(4))*t0y(nn(2),nn(5))*t0z(nn(3),nn(6)); 
       d2=t0x(nn(1),nn(4))*t2y(nn(2),nn(5))*t0z(nn(3),nn(6)); 
       d3=t0x(nn(1),nn(4))*t0y(nn(2),nn(5))*t2z(nn(3),nn(6)); 
       st(i,j)=d1+d2+d3; 
    end 
 end 
 

aq32 (Бутеноп edition)

 
 function[sstt]=aq32(num_pnts_xyz) 
  
 razm=prod(num_pnts_xyz); 
 pp=zeros(razm,3); 
 xxx=1:num_pnts_xyz(1); 
 yyy=1:num_pnts_xyz(2); 
 zzz=1:num_pnts_xyz(3); 
  
 for i=1:(razm/num_pnts_xyz(2)) 
    pp(num_pnts_xyz(2)*(i-1)+1:num_pnts_xyz(2)*i,2)=yyy; 
 end 
  
 for i=1:num_pnts_xyz(3) 
    pp(num_pnts_xyz(1)*num_pnts_xyz(2)*(i-1)+1:num_pnts_xyz(1)*num_pnts_xyz(2)*i,3)=zzz(i); 
 end 
  
 for i=1:num_pnts_xyz(3) 
    for j=1:num_pnts_xyz(1) 
        pp(num_pnts_xyz(2)*(j-1)+1+(i-1)*num_pnts_xyz(1)*num_pnts_xyz(2):num_pnts_xyz(2)*j+(i-1)*num_pnts_xyz(1)*num_pnts_xyz(2),1)=xxx(j); 
    end 
 end 
  
 t0x=sparse(q0(num_pnts_xyz(1))); 
 t0y=sparse(q0(num_pnts_xyz(2))); 
 t0z=sparse(q0(num_pnts_xyz(3))); 
 %Посчитаны интегралы сплайнов 
  
 %t1x=q1(n1);t1y=q1(n2);t1z=q1(n3); 
 t2x=sparse(q2(num_pnts_xyz(1))); 
 t2y=sparse(q2(num_pnts_xyz(2))); 
 t2z=sparse(q2(num_pnts_xyz(3))); 
 %Посчитаны вторые производные сплайнов 
  
 sstt=t2x(pp(:,1),pp(:,1)).*t0y(pp(:,2),pp(:,2)).*t0z(pp(:,3),pp(:,3)); 
 sstt=sstt+t0x(pp(:,1),pp(:,1)).*t2y(pp(:,2),pp(:,2)).*t0z(pp(:,3),pp(:,3)); 
 sstt=sstt+t0x(pp(:,1),pp(:,1)).*t0y(pp(:,2),pp(:,2)).*t2z(pp(:,3),pp(:,3)); 
 

 Отправлено: 31.05.10 22:28. Заголовок: Да, кстати, в моей в..

Да, кстати, в моей версии количество узлов по x и по у могут быть разными.

 Отправлено: 21.06.10 12:54. Заголовок: А ты пригласишь нас ..

А ты пригласишь нас на своё первое преподавание?

 Отправлено: 26.12.10 00:07. Заголовок: А вот и последняя ве..

А вот и последняя версия, в связи с познанием функции kron

function[sstt]=aq32_beta_kron(num_pnts_xyz)

t0x=sparse(q0(num_pnts_xyz(1)));
t0y=sparse(q0(num_pnts_xyz(2)));
t0z=sparse(q0(num_pnts_xyz(3)));

t2x=sparse(q2(num_pnts_xyz(1)));
t2y=sparse(q2(num_pnts_xyz(2)));
t2z=sparse(q2(num_pnts_xyz(3)));

tmp=kron(t2x,t0y);
sstt=kron(t0z,tmp);

tmp=kron(t0x,t2y);
sstt=sstt+kron(t0z,tmp);

tmp=kron(t0x,t0y);
sstt=sstt+kron(t2z,tmp);

 Отправлено: 27.04.11 15:02. Заголовок: Александр, браво!..

Александр, браво!