Форум » Программистские штучки » Стабилизатор в R3 » Ответить
Стабилизатор в R3
Прохожий: По многочисленным просьбам трудящихся ниже приведён код двух вариантов счёта стабилизирующего функционала в пространстве R3. Первый-исходный, написанный Андреем Михайловичем Волковым. Второй мой. Суть у них абсолютно одинаковая, но реализация существенно отличается. aq32.m (Волков edition) [pre2] 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 [/pre2] aq32 (Бутеноп edition) [pre2] 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)); [/pre2]
Ответов - 4
Прохожий: Да, кстати, в моей версии количество узлов по x и по у могут быть разными.
D: А ты пригласишь нас на своё первое преподавание?
Прохожий: А вот и последняя версия, в связи с познанием функции 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);
administrator: Александр, браво!
полная версия страницы