Дополнительные примеры
Матричные комплекснозначные вычисления
Создаем случайную комплекснозначную матрицу и вычисляем несколько значений функции Бесселя.
>> x=rand(5,5)+(0+1i)*rand(5,5)
x =
0.6649 + 0.6020i 0.6739 + 0.4222i 0.5485 + 0.8580i 0.7009 + 0.4983i 0.6343 + 0.8983i
0.3654 + 0.2536i 0.9994 + 0.9614i 0.2618 + 0.3358i 0.9623 + 0.4344i 0.8030 + 0.7546i
0.1400 + 0.8735i 0.9616 + 0.0721i 0.5973 + 0.6802i 0.7505 + 0.5625i 0.0839 + 0.7911i
0.5668 + 0.5134i 0.0589 + 0.5534i 0.0493 + 0.0534i 0.7400 + 0.6166i 0.9455 + 0.8150i
0.8230 + 0.7327i 0.3603 + 0.2920i 0.5711 + 0.3567i 0.4319 + 0.1133i 0.9159 + 0.6700i
>> bessel(0.6,x)
ans =
0.6513 + 0.2129i 0.5987 + 0.1500i 0.7188 + 0.3544i 0.6288 + 0.1674i 0.7707 + 0.3259i
0.4243 + 0.1477i 0.8920 + 0.1439i 0.3871 + 0.2205i 0.6758 + 0.0800i 0.7521 + 0.2024i
0.5298 + 0.5571i 0.6231 + 0.0137i 0.6576 + 0.2658i 0.6641 + 0.1711i 0.4623 + 0.5360i
0.5860 + 0.2142i 0.3515 + 0.4157i 0.1348 + 0.0727i 0.6799 + 0.1904i 0.8091 + 0.1503i
0.7475 + 0.1884i 0.4303 + 0.1697i 0.5443 + 0.1510i 0.4387 + 0.0612i 0.7420 + 0.1374i
Немного меняем ее минор поэлементным умножением:
>> x(2:3,2:3)=[1 2;3 4];x(2:3,2:3)=x(2:3,2:3).*x(2:3,2:3)
x =
0.6649 + 0.6020i 0.6739 + 0.4222i 0.5485 + 0.8580i 0.7009 + 0.4983i 0.6343 + 0.8983i
0.3654 + 0.2536i 1.0000 4.0000 0.9623 + 0.4344i 0.8030 + 0.7546i
0.1400 + 0.8735i 9.0000 16.0000 0.7505 + 0.5625i 0.0839 + 0.7911i
0.5668 + 0.5134i 0.0589 + 0.5534i 0.0493 + 0.0534i 0.7400 + 0.6166i 0.9455 + 0.8150i
0.8230 + 0.7327i 0.3603 + 0.2920i 0.5711 + 0.3567i 0.4319 + 0.1133i 0.9159 + 0.6700i
Решения
Практическое занятие №2, п.1.
3.
% Книга расположена в текущем каталоге 21.xls
%Для корректного задания ссылок проверить в Excel(2003) стиль ссылок
%Убрать флаг Excel->Сервис->Параметры->Общие->Стиль ссылок R1C1
>> W=xlsread('21.xls', 1, 'A1:I5')
W =
1 1 2 NaN 1 3 -1 0 1
3 5 5 NaN 1 -5 2 3 4
-1 2 -2 NaN 2 5 -2 1 11
0 3 1 NaN NaN NaN NaN NaN NaN
1 4 -11 NaN NaN NaN NaN NaN NaN
>> B=W(:,1:3),A=W(1:3,5:9)
B =
1 1 2
3 5 5
-1 2 -2
0 3 1
1 4 -11
A =
1 3 -1 0 1
1 -5 2 3 4
2 5 -2 1 11
4.
>> A(2,3)=x*A(2,3)-1, B(3,1)=sin(y)*B(3,1)+1
5.
>> A.*B'
ans =
1.0000 9.0000 -0.3530 0 1.0000
1.0000 -25.0000 -5.4000 9.0000 16.0000
4.0000 25.0000 4.0000 1.0000 -121.0000
6.
>> v=linspace(1,1,2);C=diag([v 1])+diag(x*v,1)+diag(y*v,-1),
C =
1.0000 -0.8500 0
9.7856 1.0000 -0.8500
0 9.7856 1.0000
>> A*B+C
ans =
11.6470 17.1500 8.0000
-1.1676 -3.4000 -59.4500
27.2939 79.7856 -86.0000
7.
>> dlmwrite('myfile.txt', ans, 'delimiter', '\t', 'precision', 6)
8.
>> trace(inv(ans))
ans =
0.1305
Практическое занятие №2, п.2.
1.
% Ваши численные результаты в лабораторной работе могут отличаться
% ввиду случайности значений элементов матриц. Для самопроверки можно
% принудительно ввести именно данные матрицы А и В
>> A=rand(3)+i*rand(3), B=expm(A), z=B\ones(3,1)
A =
0.9501 + 0.4447i 0.4860 + 0.9218i 0.4565 + 0.4057i
0.2311 + 0.6154i 0.8913 + 0.7382i 0.0185 + 0.9355i
0.6068 + 0.7919i 0.7621 + 0.1763i 0.8214 + 0.9169i
B =
0.7557 + 1.3942i -1.0258 + 2.5609i -1.2289 + 0.8812i
-1.7591 + 1.0176i 0.0147 + 1.6752i -2.2026 + 0.9689i
-1.0854 + 2.3869i -0.8267 + 1.9390i -0.1626 + 2.1086i
z =
-0.0997 - 0.1304i
0.0039 - 0.3194i
-0.0813 + 0.0232i
2.
% Тех же результатов можно было бы достичь через команду diag или даже без нее
>> C=cat(1,cat(2,A*B,B),cat(2,A,A.*B));D1=real(C);D2=imag(C)
D1 =
-3.1588 -4.8146 -4.4527 0.7557 -1.0258 -1.2289
-5.2554 -4.8658 -5.4803 -1.7591 0.0147 -2.2026
-5.2456 -5.3915 -5.3599 -1.0854 -0.8267 -0.1626
0.9501 0.4860 0.4565 0.0980 -2.8592 -0.9185
0.2311 0.8913 0.0185 -1.0329 -1.2235 -0.9472
0.6068 0.7621 0.8214 -2.5489 -0.9718 -2.0669
D2 =
1.1829 3.3544 -0.3721 1.3942 2.5609 0.8812
-0.5754 0.7272 -1.4280 1.0176 1.6752 0.9689
2.8754 2.8557 1.4947 2.3869 1.9390 2.1086
0.4447 0.9218 0.4057 1.6607 0.2990 -0.0963
0.6154 0.7382 0.9355 -0.8474 1.5040 -2.0425
0.7919 0.1763 0.9169 0.5889 1.3320 1.5829
3.
>> Electr=@(x,y) sum(sum(((x-D1).^2+(y-D2).^2).^(-0.5))); Electr(0,0)
ans =
21.4684
4.
% при нехватке памяти/мощности подождать, уменьшить число точек до 100 или 500
>> [X,Y]=meshgrid(linspace(-2,2,1000));
>> Z=[];for k=1:length(X), for n=1:length(Y), Z(k,n)=Electr(X(1,k),Y(n,1)); end ; end
>> contourf(X,Y,Z, logspace(0,1.8,10));colorbar;
5.
>> h=4/1000;[U,V] = gradient(Z,h,h);
% Делаем матрицы разреженными, чтобы стрелки не сливались
>> N=20;NN=1000/N;x=[];y=[];z=[];u=[];v=[];
>> for i=1:N, for j=1:N, m=NN*(i-1)+1;n=NN*(j-1)+1;…
x(i,j)=X(m,n);y(i,j)=Y(m,n);z(i,j)=Z(m,n);u(i,j)=U(m,n);v(i,j)=V(m,n); end; end
% Строим лишний график ради осей полярных координат
h = polar([0 2*pi], [0 2]); delete(h);hold on;
% График скоростей и контурный, т.е. напряженностей и потенциала
h=quiver(gca,x,y,z,v,8.0);set(h,'Color','r');contour(x,y,z);legend show;
6.
% Строим логарифм потенциала, а цвет рациональнее изменять по обратной величине
>>surf(X,Y,log10(Z),Z.^(-1),'LineStyle','none');colormap hsv;colorbar
7.
% Оформим в виде М-файла.
% Обратите внимание, что функция Electr вводится внутри цикла
% Потенциал нормируется на минимальное значение, и строится его логарифм
function M=MyCin(a)
tt=0:0.05:2*pi;[x,y]=meshgrid(linspace(-2,2,501));z=[];MM=[];
a=0.5;I=ones(6,1)*linspace(1,6,6); J=I';
D1=6*rand(6)-3;D2=6*rand(6)-3; D10=D1;D20=D2;
flex1=@(t) D10.*(1+a*cos(I+t*J));
flex2=@(t) D20.*(1+a*sin(I+t*J));
Hfig=figure(1);Hax=axes();Hsur=surf(x,y,ones(length(x)),'LineStyle','none');
axis manual;axis([-2 2 -2 2 0 6]);caxis([0 6]);title(''); colorbar; colormap jet;
for t=tt,
D1=flex1(t);D2=flex2(t);
Electr=@(x,y) sum(sum(((x-D1).^2+(y-D2).^2).^(-0.5)));
for k=1:length(x),
for n=1:length(y),
z(k,n)=Electr(x(1,k),y(n,1));
end ;
end;
z=z-min(min(z))+1;title(strcat('Время: ',num2str(t),' сек'));
set(Hsur,'ZData',log(z),'CData',log(z)); drawnow;
MM=[MM,getframe(Hfig)];
end;
movie2avi(MM,'mym1.avi','fps',4);
M=1;
end
Практическое занятие №3, п.1.
1.
function Lab3=Lab31(A)
rrr=@(t) cos(A*t).^2;
zzz=@(t) 1+sqrt(A*rrr(t)).*sin(t);
fi=linspace(0,2*pi,201);ro=rrr(fi);Z=zzz(fi);
[X,Y,Z]=pol2cart(fi,ro,Z);figure(1);plot3(X,Y,Z);grid on;
figure(2);comet3(X,Y,Z);Lab3=0;
end
 
2.
function res=myfun1(a,s)
if ((isvector(a))&(ischar(s)))
switch s(1)
case 'a'
w=sprintf('Среднее значение: %f',mean(a));
case 'p'
w=sprintf('Полином в точке: %f',polyval(a,1.1));
case 'g'
w=sprintf('График');plot(a);
otherwise
w=sprintf('Неверные данные!');
end;
disp(w);res=0;
else display('Неверный формат данных!');res=1;
end;
end
>> t=10*rand(10,1)-6;myfun1(t,'pbx');
Полином в точке: -50.535823
3.
%Подготовьте xls-файл, столбцы А,В,С – предмет, имя, оценка
% Массив ячеек D имеет ту же структуру. Шапки не делать
function M=myfun2(flag)
global Dbase
global Dip
function R=reading()
%Заполните D сами
R=0;
end
D={};
%Определение способа ввода - ручного или из Ексел-файла
if (isnumeric(flag))
reading();
else
[D1,D2]=xlsread('1234.xls',1,'A:C');D=cat(2, D2, num2cell(D1));clear 'D2' 'D1';
end
%D-столбцы данных, их в массив записей Dbase
Dbase=cell2struct(D,{'Subject','Name','Value'},2);
%Данные из массива ячеек в вектора
Subject=char(D{:,1});Name=char(D{:,2}); Value=cell2mat(D(:,3)); clear('D');
[Name, Permut]=sortrows(Name);E=Dbase;
for i=1:length(E), Dbase(i)=E(Permut(i)); end; clear 'E';
%База данных пересортирована
%Теперь вычисляем среднее, делаем структуру Dip
[UnN, NN]=unique(Name,'rows');
NN=[0;NN];dip=struct('Name',{'xxx'},'Merit',{5});Dip=[];
for i=1:size(UnN,1), dip.Name=UnN(i,:);db=0;
for j=(NN(i)+1):NN(i+1),
db=db+Dbase(j).Value;
end;
dip.Merit=db/(NN(i+1)-NN(i));Dip=[Dip;dip];
end;
%Усложняя:Структуру в вектора, их - на печать
DipCell=(struct2cell(Dip))',
DipCellName=char(DipCell{:,1});DipCellVal=cell2mat(DipCell(:,2));
for i=1:size(Dip),
disp(strcat('Студент: ',DipCellName(i,:),' Балл: ',num2str(DipCellVal(i)))),
end;
M=0;
end
Практическое занятие №3, п.2.
% Перед выполнением работы составить файл MyMath.m по шаблону п.1 Примера
1.
% Внесено изменение в алгоритм. Итерации идут группами, примерно по сто.
function M1=M1(x,p,flag)
MassOfPoint=@(x) 1;
function Region=Region(x)
if (sum(x.*x)>1) Region=0; else Region=MassOfPoint(x); end
end
N=1;NN=100*2^(p);eps=p*(1e-6);
med1=0;med=1;media=[];summa=0;
while ((abs(med1-med)>eps)&&(N*sqrt(eps)<1))
for i=1:NN, ksi=rand(1,p);summa=summa+Region(ksi); end;
media(N)=summa/NN;med=med1;med1=mean(media); summa=0;
%если захотим увидеть ход процесса, каждый десятый шаг
if (flag&&(~mod(N,10))) ['процессинг: ' num2str(med1) ...
' ' num2str(med)], end;
N=N+1;
end
%вторым аргументом идет погрешность расчета
M1=[med1*((2*x).^p) 1/sqrt(N*NN)];
end
….
case '1'
tic;Z=[];Z1=[];
for n=1:K, disp(n),res=M1(1,n,false); Z=[Z;res(1),res(2)];end;
for n=1:K, disp(n),res=M1(1,n,false); Z1=[Z1;res(1),res(2)];end;
VolSp=@(x) pi^(x/2).*(1./gamma(0.5*x+1));
E=Z(:,2)'*ones(K);E1=Z1(:,2)'*ones(K);norm(Z(:,1)-Z1(:,1)),
x=1:K;figure(1);hold on;
plot(x,Z(:,1),'-r');errorbar(x,Z(:,1),E,'-r');
plot(x,Z1(:,1),'-g');errorbar(x,Z1(:,1),E1,'-g');
fplot(VolSp,[1 K],'-m');
legend('Монте-Карло1','dy','Монте-Карло2','dy','Теория');
grid on;legend show;
M=Z;toc,
>> MyMath('1',12)
…
ans =
0.3075
E lapsed time is 222.764131 seconds.
ans =
2.0000 0.0408
3.1900 0.0151
4.1597 0.0060
4.9331 0.0025
5.2497 0.0021
5.2142 0.0022
4.6321 0.0023
4.0785 0.0017
3.3158 0.0012
2.6167 0.0016
1.8250 0.0013
1.3400 0.0009
2.
function M2=M2(a,b)
w=2*pi/(b-a);
function Fsin=Fsin(x)
Fsin=P5(x).*sin(w*n*x);
end;
function Fcos=Fcos(x)
Fcos=P5(x).*cos(w*n*x);
end;
for n=0:10,
A(n+1)=quadl(@Fsin,a,b);B(n+1)=quadl(@Fcos,a,b);
end;
B(1)=B(1)/2;M2=abs(A+i*B)*2/(b-a);
end
…
P5=@(x) polyval([5 0 4 3 2 1],x);
switch
…
case '2'
roots([5 0 4 3 2 1]),x0=[-1 -0.5 0.5 1];[x0;P5(x0)],
cf=abs(fft(P5(linspace(-1,1,11))'));tf=M2(-1,1)';
[cf tf],
M=2;
>> MyMath('2',12)
ans =
0.4421 + 0.9076i
0.4421 - 0.9076i
-0.1798 + 0.5845i
-0.1798 - 0.5845i
-0.5247
ans =
-1.0000 -0.5000 0.5000 1.0000
-7.0000 0.0938 3.4063 15.0000
ans =
24.2000 2.0000
22.9582 3.1683
21.3636 2.6417
18.1306 2.0586
15.9210 1.6313
14.9136 1.3387
14.9136 1.1311
15.9210 0.9776
18.1306 0.8601
21.3636 0.7673
22.9582 0.6924
% Расхождение левой и правой колонок обусловлено несовпадением частоты
% т.е. особенностью вычисления быстрого Фурье-преобразования
3.
%%
%Тестовая функция
function M3=M3(x,y,z)
A=K(1,1);h=1e-4;
m3=((x.^2).*y).*z;
function Dz=Dz(x,y,z)
part=@(x,z) A*x.*(z.^2);
partZ=cat(4,part(x,z-h),part(x,z),part(x,z+h));
% Матрица по 4-му измерению, ведь x,y,z - могут быть трехмерными
% Градиент по переменной берется по 2-му измерению
partZ=permute(partZ,[1 4 3 2]);
Dz=gradient(partZ,h);Dz=permute(Dz,[1 4 3 2]);
Dz=Dz(:,:,:,2);Dz=reshape(Dz,size(m3));
end
M3=m3-(1+x.^2).*Dz(x,y,z);
end
%%
%Рабочая функция
function M4=M4(x,y,z)
A=K(2,1);B=K(2,2);h=1e-4;
m4=A*P5(x.*y)+erf(x-y+B.*z);
function Dz=Dz(x,y,z)
function part0=part0(x,y,z)
p=legendre(5,sin(x.*z+A*y));p=p(1,:,:,:);
part0=reshape(p,size(m4));
end;
part=@(x,y,z) z.*part0(x,y,z);
partZ=cat(4,part(x,y,z-h),part(x,y,z),part(x,y,z+h));
partZ=permute(partZ,[1 4 3 2]);
Dz=gradient(partZ,h);Dz=permute(Dz,[1 4 3 2]);
Dz=reshape(Dz(:,:,:,2),size(m4));
end;
M4=m4-B*Dz(x,y,z);
end
%%
%Строим графики по заданному указателю функции ММ
function M5=M5(MM,Lim,Name,flg)
%Ниже - процедура 3Д-графики, flg, если режим теста
function Mygraph3d=Mygraph3d(X,Y,Z,U,aver,s1,flag)
if flag
figure('Name',['3D-графики: ' Name]);axis(Lim);grid on;box on;
xlabel('x');ylabel('y');zlabel('z');
view(3);light('Position',[Lim(2) Lim(4) Lim(6)]);camlight;lighting flat;
else
hold on;
end;
p = patch(isosurface(X,Y,Z,U,aver));isonormals(X,Y,Z,U,p);
set(p,'FaceColor',s1,'EdgeColor','none');
Mygraph3d=gca;
end;
%Построение линий уровня на 2D
function Mygraph2d=Mygraph2d(level)
%level - либо задаем число линий, либо вектор значений уровней
Hfig=figure('Name',['Семейство линий уровня для лаплассиана ' Name]);
axis([Lim(1) Lim(2) Lim(5) Lim(6)]);hold on;
xlabel('x');ylabel('z');LSpec={'-r',':b','-.g'};
%Фиксируем y, строим линии уровня LU(x,const,z)=const
[X,Z]=meshgrid(x,z);
for yy=1:3,
n=1+ceil((yy-1)*N/2);Ly=squeeze(L(n,:,:))';
[C,Hcg]=contour(gca,X,Z,Ly,level,LSpec{yy});
set(Hcg,'DisplayName',['at y=' num2str(y(n))],'ShowText','off');
if (yy==1) set(Hcg,'ShowText','on'); end;
end;
grid on;legend show;
Mygraph2d=Hfig;
end;
%Расчеты
N=100;x=linspace(Lim(1),Lim(2),N+1);y=linspace(Lim(3),Lim(4),N+1);
z=linspace(Lim(5),Lim(6),N+1);[X,Y,Z]=meshgrid(x,y,z);U=MM(X,Y,Z);
L=6*del2(U,(Lim(2)-Lim(1))/N,(Lim(4)-Lim(3))/N,(Lim(6)-Lim(5))/N);
Um=mean(mean(mean(U)));Lm=mean(mean(mean(L)));
s{1}=['U(r)=',num2str(Um)];s{2}=['d2U(r)=',num2str(Lm)];
%Данные и лаплассиан сформированы. Далее графика
if flg
Mygraph3d(X,Y,Z,U,Um,'red',true);title('Тестовый пример - функция и ее лаплассиан');
Mygraph3d(X,Y,Z,L,Lm,'blue',false);
legend show;legend(s,'Location','NorthEast');
hold off;
else
Mygraph3d(X,Y,Z,U,Um,'green',true);title('Рабочая функция');
legend show;legend(s{1},'Location','NorthEast');
Mygraph3d(X,Y,Z,L,Lm,'magenta',true);title('Лаплассиан рабочей функции');
legend show;legend(s{2},'Location','NorthEast');
end;
%Построены две поверхности. Теперь строим линии уровня лаплассиана
Mygraph2d(3);
M5=0;
end
%%
%Тестовая функция
function M3=M3(x,y,z)
A=K(1,1);h=1e-4;
m3=((x.^2).*y).*z;
function Dz=Dz(x,y,z)
part=@(x,z) A*x.*(z.^2);
partZ=cat(4,part(x,z-h),part(x,z),part(x,z+h));
% Матрица по 4-му измерению, ведь x,y,z - могут быть трехмерными
% Градиент по переменной берется по 2-му измерению
partZ=permute(partZ,[1 4 3 2]);
Dz=gradient(partZ,h);Dz=permute(Dz,[1 4 3 2]);
Dz=Dz(:,:,:,2);Dz=reshape(Dz,size(m3));
end
M3=m3-(1+x.^2).*Dz(x,y,z);
%M3=A*sin(z)-x.^3-y.^4;
end
%%
%Рабочая функция
function M4=M4(x,y,z)
A=K(2,1);B=K(2,2);h=1e-4;
m4=A*P5(x.*y)+erf(x-y+B.*z);
function Dz=Dz(x,y,z)
function part0=part0(x,y,z)
p=legendre(5,sin(x.*z+A*y));p=p(1,:,:,:);
part0=reshape(p,size(m4));
end;
part=@(x,y,z) z.*part0(x,y,z);
partZ=cat(4,part(x,y,z-h),part(x,y,z),part(x,y,z+h));
partZ=permute(partZ,[1 4 3 2]);
Dz=gradient(partZ,h);Dz=permute(Dz,[1 4 3 2]);
Dz=reshape(Dz(:,:,:,2),size(m4));
end;
M4=m4-B*Dz(x,y,z);
%M4=P5(x.*y);
end
%%
%Строим графики по заданному указателю функции ММ
function M5=M5(MM,Lim,Name,flg)
%Ниже - процедура 3Д-графики, flg, если режим теста
function Mygraph3d=Mygraph3d(X,Y,Z,U,aver,s1,flag)
if flag
figure('Name',['3D-графики: ' Name]);axis(Lim);grid on;box on;
xlabel('x');ylabel('y');zlabel('z');
view(3);light('Position',[Lim(2) Lim(4) Lim(6)]);camlight;lighting flat;
else
hold on;
end;
p = patch(isosurface(X,Y,Z,U,aver));isonormals(X,Y,Z,U,p);
set(p,'FaceColor',s1,'EdgeColor','none');
Mygraph3d=gca;
end;
%Построение линий уровня на 2D
function Mygraph2d=Mygraph2d(level)
%level - либо задаем число линий, либо вектор значений уровней
Hfig=figure('Name',['Семейство линий уровня для лаплассиана ' Name]);
axis([Lim(1) Lim(2) Lim(5) Lim(6)]);hold on;
xlabel('x');ylabel('z');LSpec={'-r',':b','-.g'};
%Фиксируем y, строим линии уровня LU(x,const,z)=const
[X,Z]=meshgrid(x,z);
for yy=1:3,
n=1+ceil((yy-1)*N/2);Ly=squeeze(L(n,:,:))';
[C,Hcg]=contour(gca,X,Z,Ly,level,LSpec{yy});
set(Hcg,'DisplayName',['at y=' num2str(y(n))],'ShowText','off');
if (yy==1) set(Hcg,'ShowText','on'); end;
end;
grid on;legend show;
Mygraph2d=Hfig;
end;
%Расчеты
N=100;x=linspace(Lim(1),Lim(2),N+1);y=linspace(Lim(3),Lim(4),N+1);
z=linspace(Lim(5),Lim(6),N+1);[X,Y,Z]=meshgrid(x,y,z);U=MM(X,Y,Z);
L=6*del2(U,(Lim(2)-Lim(1))/N,(Lim(4)-Lim(3))/N,(Lim(6)-Lim(5))/N);
Um=mean(mean(mean(U)));Lm=mean(mean(mean(L)));
s{1}=['U(r)=',num2str(Um)];s{2}=['d2U(r)=',num2str(Lm)];
%Данные и лаплассиан сформированы. Далее графика
if flg
Mygraph3d(X,Y,Z,U,Um,'red',true);title('Тестовый пример - функция и ее лаплассиан');
Mygraph3d(X,Y,Z,L,Lm,'blue',false);
legend show;legend(s,'Location','NorthEast');
hold off;
else
Mygraph3d(X,Y,Z,U,Um,'green',true);title('Рабочая функция');
legend show;legend(s{1},'Location','NorthEast');
Mygraph3d(X,Y,Z,L,Lm,'magenta',true);title('Лаплассиан рабочей функции');
legend show;legend(s{2},'Location','NorthEast');
end;
%Построены две поверхности. Теперь строим линии уровня лаплассиана
Mygraph2d(3);
M5=0;
end
…
case '3'
tic;v=[1 2 5 10 0 2];M5(@M3,v,'Test',true);
v=[-1 1 0 1 -1 0];M5(@M4,v,'Work',false);
M=3;toc,
>> K=[1 -1;0.1 -1];
>> MyMath('3',K);
Elapsed time is 30.782295 seconds.
4.
%%
function M6=M6(MM)
N=101;scrsz = [1 1 0.8 0.8].*get(0,'ScreenSize');s='';
figure('Name','Simulation Plot Window','NumberTitle','off','Position',scrsz);
%1
subplot(2,2,1);x=linspace(K(3,1)/2,2*K(3,1),N);
[s,error]=sprintf('U(x,y0,z0) at y0=%3.2f, z0=%3.2f',K(3,2),K(3,3));
if (func2str(MM)=='MyMath/M3')
[leg,error]=sprintf('a=%3.2f',K(1,1));
else
[leg,error]=sprintf('(a,b)=(%3.2f,%3.2f)',K(2,1),K(2,2));
end
U=MM(x,K(3,2),K(3,3));plot(x,U,'r','DisplayName',leg);
title(s);grid on;legend('Location','NE');legend show;xlabel('x');ylabel('U(x)');
%2
subplot(2,2,2);y=linspace(K(3,2)/2,2*K(3,2),N);[X,Y]=meshgrid(x,y);U=MM(X,Y,K(3,3));
[s,error]=sprintf('U(x,y,z0) at z0=%3.2f',K(3,3)); contourf(x,y,U,'DisplayName',leg);
title(s);xlabel('x');ylabel('y');grid on;legend('Location','NE');legend show;colorbar;
%3
subplot(2,2,3);z=linspace(K(3,3)/2,2*K(3,3),N);[X,Y,Z]=meshgrid(x,y,z);U=MM(X,Y,Z);
d=linspace(0.5,2,4);xd=K(3,1)*d;yd=K(3,2)*d;zd=K(3,3)*d(2:3);
h = slice(x,y,z,U,xd,yd,zd);set(h,'FaceColor','interp','EdgeColor','none');
axis tight;title('Объемный график сечений для U(x,y,z)');
xlabel('x');ylabel('y');zlabel('z');grid on;colorbar;
%4
PU=[];
if (func2str(MM)=='MyMath/M3')
[s,error]=sprintf('f(x,a)=U(x,y0,z0) at y0=%3.2f, z0=%3.2f',K(3,2),K(3,3));
a=K(1,1);y=K(3,2);z=K(3,3);aa=linspace(a/2,2*a,N);
for k=aa, K(1,1)=k;d=MM(x,y,z);PU=[PU;d]; end;K(1,1)=a;
[X,Y]=meshgrid(x,aa);subplot(2,2,4);surfc(X,Y,PU,'LineStyle','None');
xlabel('x');ylabel('a');zlabel('U(x,a)');colorbar;axis tight;axis on;box on;
title(s);
else
[s,error]=sprintf('f(a,b)=U(x0,y0,z0) at x0=%3.2f, y0=%3.2f, z0=%3.2f',K(3,1),K(3,2),K(3,3));
a=K(2,1);b=K(2,2);x=K(3,1);y=K(3,2);z=K(3,3);
aa=linspace(a/2,2*a,N);bb=linspace(b/2,2*b,N);
for m=aa, d=[];K(2,1)=m;for n=bb, K(2,2)=n;d=[d;MM(x,y,z)]; end; PU=[PU d]; end;
K(2,1)=a;K(2,2)=b;
[X,Y]=meshgrid(aa,bb);subplot(2,2,4);surfc(X,Y,PU,'LineStyle','None');
xlabel('a');ylabel('b');zlabel('U(a,b)');colorbar;axis tight;axis on;box on;
title(s);
end
M6=0;
end
…
case '4'
tic;M6(@M3);M6(@M4);M=4;toc,
>> K=[2 0 0;0.2 1 0;-0.5 0.5 1];
>> MyMath('4',K);
Elapsed time is 36.239721 seconds.
5.
case '5'
N=101;x=linspace(K(2,1),K(2,2),N);U=M3(x,K(3,2),K(3,3))-0.1*P5(x);
plot(x,U);grid on;
%повторить несколько раз предыдущие две строки
M=5;
>> K=[2 0 0;-1 2 0;0 0.5 1]
K =
2.0000 0 0
-1.0000 2.0000 0
0 0.5000 1.0000
>> MyMath('5',K);
>> K(2,1)=-0.1;K(2,2)=0.1;
>> MyMath('5',K);
>> K(2,1)=-0.03;K(2,2)=-0.02;
>> MyMath('5',K);
>> K(2,1)=-0.024;K(2,2)=-0.023;
>> K(2,1)=-0.024;K(2,2)=-0.023;
>>
>> MyMath('5',K);
>> K(2,1)=-0.0238;K(2,2)=-0.0236;
>> MyMath('5',K);
>> x=-0.0238;
…
U=M3(x,K(3,2),K(3,3));q=polyfit(x,U,5),q=q-0.1*[5 0 4 3 2 1];roots(q),
M=5;
>> K=[2 0 0;-1 2 0;0 0.5 1]
K =
2.0000 0 0
-1.0000 2.0000 0
0 0.5000 1.0000
>> MyMath('5',K);
q =
-0.0000 0.0000 -4.0000 0.5000 -4.0000 0.0000
ans =
-0.0321 + 2.7775i
-0.0321 - 2.7775i
0.0440 + 1.0434i
0.0440 - 1.0434i
-0.0238
% Нетрудно видеть, что корни и коэффициенты полинома найдены верно
Практическое занятие №4
2.
Для масштабируемости окна и его элементов устанавливаем GUI-options значение свойства масштабируемости Resizable, но сверх того во всех элементах в GUIDE свойство Units как Proportional. Средствами GUIDE размечаем саму figure, оси, статический текст, текст в боксе и кнопку пуска. Средствами GUIDE размечаем само меню, затем сохраняем оба файла (lab42.fig и lab42.m). Далее редактируем последний файл самостоятельно, чтобы: во-первых, доопределить обработчики, а во-вторых, создать логотип, панель инструментов и задать их обработчики.
Листинг файла lab42.m без избыточных комментариев и процедур см. ниже:
function varargout = lab42(varargin)
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @lab42_OpeningFcn, ...
'gui_OutputFcn', @lab42_OutputFcn, ...
'gui_LayoutFcn', [], ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before lab42 is made visible.
function lab42_OpeningFcn(hObject, eventdata, handles, varargin)
% Choose default command line output for lab42
handles.output = hObject;
%Добавление элементов
handles.toolbar1=uitoolbar('Tag','toolbar1');
handles.tbr_button1=uipushtool(handles.toolbar1,'Tag','tbr_button1');
handles.tbr_button2=uipushtool(handles.toolbar1,'Tag','tbr_button2');
handles.tbr_button3=uipushtool(handles.toolbar1,'Tag','tbr_button3');
%1-й способ - создание картинки на лету
x=ones(20);CData=cat(3,x,x/2,x/3);set(handles.tbr_button1,'CData',CData);
%2-й способ - чтение из bmp-файла. Функция imresize не из ядра МАТЛАБ
x=imread('photo.bmp','bmp');CData=imresize(x,[20 20],'bilinear');
set(handles.tbr_button2,'CData',CData,'TooltipString','2-я кнопка');
%3-й способ - чтение из ico-файла с помощью утилиты iconRead
% работа утилиты зависит от формата файла и не всегда корректна
icopath=fullfile(matlabroot,'help','techdoc','creating_guis','examples');
path(icopath,path);CData=iconRead('my.ico');rmpath(icopath);
set(handles.tbr_button3,'CData',CData,'TooltipString','3-я кнопка');
%Поместим логотип в верхнем правом углу
mn=[50 50];x=imread('photo.bmp','bmp');CData=imresize(x,mn,'bilinear');
h=handles.output;set(h,'Units','pixels');pos=get(h,'Position');
pos=[pos(3)-mn(1) pos(4)-mn(2) mn];
haxes_im=axes('Tag','axes_im','Units','pixels','Position',pos,'Parent',h);
set(h,'Units','normalized');set(haxes_im,'Units','normalized');
handles.haxes_im=haxes_im;handles.im=image(CData,'Parent',haxes_im,'Tag','im');
set(haxes_im,'Visible','off');
%Установка вызова обработчика для первой кнопки панели инструментов
set(handles.tbr_button1,'ClickedCallback',{@tbr1_Callback, handles,'Апельсины'});
% Update handles structure
guidata(hObject, handles);
% --- Outputs from this function are returned to the command line.
function varargout = lab42_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;
% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
axes(handles.axes1);cla reset;set(gca,'Color','Yellow');
text(0.1,0.1,'Да здравствует Лимон!');
% --------------------------------------------------------------------
function Kapusta_Callback(hObject, eventdata, handles)
% hObject handle to Kapusta (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
beep;
function tbr1_Callback(hObject, eventdata, handles, label)
axes(handles.axes1);cla reset;set(gca,'Color',[1 1/2 1/3]);
text(0.1,0.1,label);title(label);set(hObject,'ToolTipString','Apelsin');
Список литературы
Дьяконов В. П. Справочник по применению системы PC MATLAB. — М.: «Физматлит», 1993. — С. 112.
Иглин С.П. Математические расчёты на базе MATLAB. БХВ, 2005, Санкт-Петербург, Россия, 640 с.
Дьяконов В.П., Круглов В. Математические пакеты расширения MATLAB. Специальный справочник. – Спб.: Питер, 2001. – 480 с.
Дьяконов В.П., Круглов В. MATLAB. Анализ, идентификация и моделирование систем. Специальный справочник. – Спб.: Питер, 2002. – 448 с.
Дэбни Дж., Хартман Т. Simulink 4. Секреты мастерства. – М.: БИНОМ. Лаборатория знаний, 2003. – 403 с.
Загидуллин Р.Ш. LabView в исследованиях и разработках. – М.: Горячая линия – Телеком, 2005. – 352 с.
Оглавление
Введение 3
Лабораторное занятие №1-№5 4
Знакомство с пакетом MATLAB 7.2. 4
1. Структура пакета и принципы работы 4
Пример: 5
Задание: 6
Практическое занятие №6-№10 7
Проведение вычислений без М-файлов 7
1. Элементарные матричные вычисления. 7
Пример: 7
Задание 9
2. Элементарные функциональные вычисления и построение графиков. 9
Пример: 12
Задание: 17
Лабораторное занятие №11- №15 18
Проведение вычислений с помощью М-файлов 18
1. Создание и отладка М-файлов. 18
Пример: 20
Задание: 25
2. Решение некоторых стандартных математических задач 25
Пример: 26
Задание: 32
Лабораторное занятие №17-№22 34
Реализация вычислений с помощью графического интерфейса (GUI) 34
Пример: 37
Задание: 46
Дополнительные примеры 47
Матричные комплекснозначные вычисления 47
Решения 47
Практическое занятие №2, п.1. 47
3. 47
4. 48
5. 48
6. 48
7. 48
8. 48
Практическое занятие №2, п.2. 48
1. 48
2. 49
3. 49
4. 49
5. 50
6. 51
7. 51
Практическое занятие №3, п.1. 52
1. 52
2. 52
3. 52
Практическое занятие №3, п.2. 53
1. 53
2. 54
3. 55
4. 61
5. 63
Практическое занятие №4 64
2. 64
Учебное издание
Математические основы моделирования
Составитель КРЫЛОВ Александр Олегович
Редактор С.И. К о с т е р и н а
Технический редактор Г.Н. Ш а н ь к о в а
Подп. в печать
Формат 60х84 1/16. Бумага офсетная.
Печать офсетная.
Усл. п. л. . Усл. кр.-отт. .
Уч.-изд. л. .
Тираж 50 экз. Заказ С. -
__________________________________________________________________________________________________
Государственное образовательное учреждение
высшего профессионального образования
«Самарский государственный технический университет»
443100 г. Самара, ул. Молодогвардейская, 244. Главный корпус. |
