Руководство для мальчиков по Matlab и Scilab

Matlab — это язык программирования, который существует уже некоторое время и остается популярным по сей день. Scilab является братом Matlab с открытым исходным кодом, поэтому мы будем использовать оба взаимозаменяемо для наших целей. Matlab и Scilab — отличные языки для вычислительного программирования для компьютерных ученых и других специалистов, которые могут использовать возможности этого языка для быстрой обработки больших наборов данных. Matlab или Scilab должны быть широко известны благодаря всем процессам, поддающимся количественной оценке, поскольку они позволяют автоматизировать и сокращать время, делая нашу жизнь немного проще. В этой статье мы рассмотрим многочисленные примеры того, как Matlab и Scilab можно применять для решения уникальных задач.

Простые вычисления с Matlab

% Name: Roberto Carlos Baldizon Diaz
% Matlab for array calculations
A=[1, 3, 5];
B=[-3, -2, 4];
dot(A, B)
ans =
11
 
a=[0, -1, -4, -8];
b=[4, -2, -3, 24];
dot(a, b)
ans =
-178
% Matlab for array calculations with scientific application
D=[1.2; 7.8; 2.7];
V=[700; 200; 300];
disp(‘Total Mass of Components in Grams:’)
M=dot(D, V)
Total Mass of Components in Grams:
M =
3210
% Matlab for simple array algebra
 
A2=[12, 4; 3, -5];
B2=[2, 12; 0, 0];
A2*B2
ans =
24 144
 6 36
 
a2=[1, 3, 5; 2, 4, 6]’;
b2=[-2, 4; 3, 8; 12, -2]’;
a2*b2
ans =
6 19 8
 10 41 28
 14 63 48
% Matlab to find array's determinants and inverses
a3=[2, -1; 4, 5];
b3=[4, 2; 2, 1];
c3=[2, 0, 0; 1, 2, 2; 5, -4, 0];
disp(‘Determinant and Inverse(if exists) of A’)
det(a3)
inv(a3)
disp(‘Determinant and Inverse(if exists) of B’)
det(b3)
disp(‘Determinant and Inverse(if exists) of B’)
det(c3)
inv(c3)
Determinant and Inverse(if exists) of A
ans =
14
ans =
0.3571 0.0714
 -0.2857 0.1429
Determinant and Inverse(if exists) of B
ans =
0
Determinant and Inverse(if exists) of B
ans =
16.0000
ans =
0.5000 0 0
 0.6250 0 -0.2500
 -0.8750 0.5000 0.2500
clear, clc

Вычисление гидродинамики с помощью Scilab

clear
//INPUT
frul=5.5; //w/m²C
ta=118; //C
ti=30; //C
//CALCULATIONS
qu=-frul*(ti-ta);
//RESULTS
qu

Солнечные инженерные расчеты с помощью Scilab

clear
//CONSTANTS
eg=0.88; 
//Glass emittence given
b=45;
//Collector tilt in degrees assumed
s=5.67*10^-8;
//Stefan-Boltzmann Constant in W/m²K^-4
//INPUT
hw=20; 
//Wind heat transfer coeff. in W/m²C
ep=0.1; 
//Plate emittence
tpc=100; 
//Plate temperature in C
tac=40; 
//Ambient temperature in C
n=1; 
//Number of covers
//CALCULATIONS
tpm=tpc+273.15; 
//Mean plate temperature in K
ta=tac+273.15;
//Plate temperature in K
f=(1+0.089*hw-0.1166*hw*ep)*(1+0.07866*n);
//Variable for Equation 6.4.9
c=520*(1–0.000051*b²);
//Variable for Equation 6.4.9
e=0.430*(1–100/tpm);
//Variable for Equation 6.4.9
//FORMULA
ut=(n/(c/tpm*((tpm-ta)/(n+f))^e)+1/hw)^-1 + (s*(tpm+ta)*(tpm²+ta²))/(1/(ep+0.00591*n*hw)+(2*n+f-1+0.133*ep)/eg-n);
//Equation 6.4.9
//OUTPUT
ut
//Result

Дополнительные расчеты солнечной энергетики с помощью Scilab

clear
//INPUT
phi=39.7392;
n=160;
lstd=-105;
lloc=-104.9903;
timemin=630;
//CALCULATIONS
delta=23.45*sind((360/365)*(284+n));
B=360*((n-1)/365);
E_t=(229.2*.000075)+229.2*(.001868*cosd(B)..
-.032077*sind(B))-229.2*(.014615*cosd(2*B)..
+.04089*sind(2*B)); // units of min
t_lambda=4*(lstd-lloc); // units of min
TC=E_t+t_lambda; // min
t_sol=timemin+TC-60; // min
t=t_sol/60 // hr
//RESULTS
omega=(t-12)*15;
d=23.45*sind(360/365*(284+n));

Еще больше расчетов солнечной энергетики с помощью Scilab

clear
//INPUT
beta = 34; // collector tilt
gamma = 5; // collector azimuth
phi= 40.44; // latitude of Pittsburgh, PA
omega= -18.81; // 11h30 solar time
delta= 14.43; // Aug 13
//CALCULATIONS
theta = acosd((sind(phi)*sind(delta)*cosd(beta))-..
(cosd(phi)*sind(delta)*sind(beta)*cosd(gamma))+..
(cosd(phi)*cosd(delta)*cosd(beta)*cosd(omega))+..
(sind(phi)*cosd(delta)*sind(beta)*cosd(gamma)*cosd(omega))+..
(cosd(delta)*sind(beta)*sind(omega)*sind(gamma)));
//RESULTS
theta
clear
//INPUT
n=131;
lstd=-90;
lloc=-87.6298;
tsolmin=900;
//CALCULATIONS
delta=23.45*sind((360/365)*(284+n));
B=360*((n-1)/365);
E_t=(229.2*.000075)+229.2*(.001868*cosd(B)..
-.032077*sind(B))-229.2*(.014615*cosd(2*B)..
+.04089*sind(2*B)); // units of min
t_lambda=4*(lstd-lloc); // units of min
TC=E_t+t_lambda; // min
timemin=tsolmin+TC-60; // min
//RESULTS
t=timemin/60 // hr

Нефтегазовая инженерия с Matlab

%%%%%%%%%Thomas Algorithm Solver Function%%%%%%%%%%%%%
function x = TDMAsolver(a,b,c,d)
n = length(b);
c(1) = c(1) / b(1);
d(1) = d(1) / b(1);
for i = 2:n-1
temp = b(i) - a(i) * c(i-1);
c(i) = c(i) / temp;
d(i) = (d(i) - a(i) * d(i-1)) / temp;
end
d(n) = (d(n) - a(n) * d(n-1))/( b(n) - a(n) * c(n-1));
x(n) = d(n);
for i = n-1:-1:1
x(i) = d(i) - c(i) * x(i + 1);
end
%%%%%%%%%%Time Steps Plotting Program%%%%%%%%%%%%%%%%%
clear;
close all;
n=100;
dx=10:20:2000;
p=ones(1,n)*1000;
p0=3000;
pl=1000;
p=[p0,p,pl];
alpha=4;
dt=1000;
for i=1:dt
p(i+1,1)=p0;
p(i+1,2+n)=pl;
for j=2:(n+1)
p(i+1,j)=(p(i,j-1)-(2-alpha)*p(i,j)+p(i,j+1))/alpha;
end;
end;
plot(dx,p(2,2:(n+1)),'-ok');
hold on;
plot(dx,p(11,2:(n+1)),'-or');
plot(dx,p(101,2:(n+1)),'-og');
plot(dx,p(1001,2:(n+1)),'-ob');
xlabel('Distance (ft)', 'Fontsize',16);
ylabel('Pressure (psia)','Fontsize',16);

a=ones(n,1);
b=ones(n,1)*(-2);
c=ones(n,1);
d=zeros(n,1);
a(1)=0;
c(n)=0;
d(1)=-p0;
d(n)=pl;
pss=TDMAsolver(a,b,c,d);
plot(dx,pss,'-om');
legend('1 Time Step','10 Time Steps','100 Time Steps', '1000 Time Steps','Steady State');
set(gca,'Fontsize',12);

%%%%%%%%%%%%%%Plot Unstable and Steady state Results%%%%%%%%%%%%
clear;
close all;
n=100;
dx=10:20:2000;
p=ones(1,n)*1000;
p0=3000;
pl=1000;
p=[p0,p,pl];
alpha=1;
dt=3;
for i=1:dt
p(i+1,1)=p0;
p(i+1,2+n)=pl;
for j=2:(n+1)
p(i+1,j)=(p(i,j-1)-(2-alpha)*p(i,j)+p(i,j+1))/alpha;
end;
end;
plot(dx,p(dt+1,2:(n+1)),'-ok');
hold on;
xlabel('Distance (ft)', 'Fontsize',16);
ylabel('Pressure (psia)','Fontsize',16);
a=ones(n,1);
b=ones(n,1)*(-2);
c=ones(n,1);
d=zeros(n,1);
a(1)=0;
c(n)=0;
d(1)=-p0;
d(n)=pl;
pss=TDMAsolver(a,b,c,d);
plot(dx,pss,'-om');
legend('3 Time Steps','Steady State');
set(gca,'Fontsize',12);

Значение, используемое в этом вопросе, равно α = 1, что составляет половину требуемой константы для стабильности, как видно на графике, результатом такого значения являются колебания в пределах отрицательных и положительных значений. На устойчивое решение не влияет альфа, поскольку она не используется при расчете результата, для которого отношение остается прежним. В этом случае использовался временной шаг Dt=3, для которого количество колебаний и результатов меньше, что позволяет сравнить поведение обеих функций. Если бы временной шаг был больше, результаты также были бы больше, а колебания были бы более выраженными, из-за чего стационарное решение выглядело бы как прямая плоская линия.

clear;
close all;
%%%%%%%%%%%%%%%Implicit Solution%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=100;
dx=10:20:2000;
p=ones(1,n)*1000;
p0=3000;
pl=1000;
alpha=4;
dt=100;
a=ones(n,1);
b=-ones(n,1)*(2+alpha);
c=ones(n,1);
d=-p*alpha;
a(1)=0;
c(n)=0;
d(1)=-p0+d(1);
d(n)=pl+d(n);
for i=1:dt
p(i+1,:)=TDMAsolver(a,b,c,d);
d=-p(i+1,:)*alpha;
d(1)=-p0+d(1);
d(n)=-pl+d(n);
end;
plot(dx,p(2,:),'-ok');
hold on;
plot(dx,p(11,:),'-om');
plot(dx,p(101,:),'-oc');
%%%%%%%%%%%%%%%Explicit Solution%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=100;
dx=10:20:2000;
p=ones(1,n)*1000;
p0=3000;
pl=1000;
p=[p0,p,pl];
alpha=4;
dt=100;
for i=1:dt
p(i+1,1)=p0;
p(i+1,2+n)=pl;
for j=2:(n+1)
p(i+1,j)=(p(i,j-1)-(2-alpha)*p(i,j)+p(i,j+1))/alpha;
end;
end;
plot(dx,p(2,2:(n+1)),'-og');
hold on;
plot(dx,p(11,2:(n+1)),'-ob');
plot(dx,p(101,2:(n+1)),'-or');
xlabel('Distance (ft)', 'Fontsize',16);
ylabel('Pressure (psia)','Fontsize',16);
legend('Implicit dt=1','Implicit dt=10','Implicit dt=100','Explicit dt=1','Explicit dt=10','Explicit dt=100');
set(gca,'Fontsize', 12);

Решения для неявной формы следуют той же схеме, что и для явной формы, в то время как временные шаги относительно малы, а постоянная альфа достигает стабильности. Вероятно, на этих временных шагах разница настолько мала, что не влияет на форму кривой, когда явная и неявная формы решения ложатся друг на друга.

clear;
close all;
%%%%%%%%%%%%%%%Implicit Solution = 0.2%%%%%%%%%%%%%%%%%%%%%%%
n=100;
dx=10:20:2000;
p=ones(1,n)*1000;
p0=3000;
pl=1000;
alpha=0.2;
dt=50;
a=ones(n,1);
b=-ones(n,1)*(2+alpha);
c=ones(n,1);
d=-p*alpha;
a(1)=0;
c(n)=0;
d(1)=-p0+d(1);
d(n)=pl+d(n);
for i=1:dt
p(i+1,:)=TDMAsolver(a,b,c,d);
d=-p(i+1,:)*alpha;
d(1)=-p0+d(1);
d(n)=-pl+d(n);
end;
plot(dx,p(dt+1,:),'-ok');
hold on;
xlabel('Distance (ft)', 'Fontsize',16);
ylabel('Pressure (psia)','Fontsize',16);
%%%%%%%%%%%%%%%Implicit Solution = 0.02%%%%%%%%%%%%%%%%%%%%%%
n=100;
dx=10:20:2000;
p=ones(1,n)*1000;
p0=3000;
pl=1000;
alpha=0.02;
dt=5;
a=ones(n,1);
b=-ones(n,1)*(2+alpha);
c=ones(n,1);
d=-p*alpha;
a(1)=0;
c(n)=0;
d(1)=-p0+d(1);
d(n)=pl+d(n);
for i=1:dt
p(i+1,:)=TDMAsolver(a,b,c,d);
d=-p(i+1,:)*alpha;
d(1)=-p0+d(1);
d(n)=-pl+d(n);
end;
plot(dx,p(dt+1,:),'-or');
%%%%%%%%%%%%%%%Explicit Solution%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=100;
dx=10:20:2000;
p=ones(1,n)*1000;
p0=3000;
pl=1000;
p=[p0,p,pl];
alpha=4;
dt=1000;
for i=1:dt
p(i+1,1)=p0;
p(i+1,2+n)=pl;
for j=2:(n+1)
p(i+1,j)=(p(i,j-1)-(2-alpha)*p(i,j)+p(i,j+1))/alpha;
end;
end;
plot(dx,p(1001,2:(n+1)),'-ob');
xlabel('Distance (ft)', 'Fontsize',16);
ylabel('Pressure (psia)','Fontsize',16);