Курсовой проект Рындин / 9206_КР_Рындин
.pdf
( ) = ln(Ф ( ) ( )) → ( ) = Ф ( ) + ( )
→ Ф ( ) = 0 → Ф ( ) = 1
Тогда, ( ) = −ln(Ф ( ) − ( )) + → Ф ( ) = ,
где U – внешнее напряжение.
Следовательно,
Ф ( ) = 1
Ф ( ) = 1
Ф ( ) = −
{ Ф ( ) =
Таким образом,
( ) = ( ( ) + √( ( ))2 + 1) 2 2
|
|
( |
) |
|
( |
) |
2 |
|
( |
) = − ( |
|
|
+ √( |
|
|
) |
+ 1) + |
|
|
|
|
|||||
|
|
2 |
|
|
2 |
|
|
|
|
|
|
|
|
|
|
||
Поиск решения системы уравнений можно осуществить с помощью метода конечных разностей, который предполагает дискретизацию дифференциальных уравнений на прямоугольных координатных сетках, то есть на сетках, элементарные ячейки которых представляют собой прямоугольники для двух измерений или параллелепипеды для трех измерений.
Таким образом, уравнение Пуассона по методу конечных разностей имеет вид:
|
2 |
|
|
+1 |
− |
|
|
|
− |
−1 |
|
|
|
|
|
|
( |
|
|
− |
|
|
) = Ф − Ф |
|
|
− − |
|||
|
|
|
|
|
|
|
|
|
|
|
||||
∆ + ∆ |
|
|
∆ |
|
|
∆ |
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
||||||
|
−1 |
|
|
|
|
|
|
|
−1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
21 |
|
|
|
|
|
Уравнения непрерывности: |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||||||
|
|
|
2 |
|
|
|
|
|
|
|
|
|
+1 + |
|
|
|
|
Ф |
|
|
− Ф |
|
|
|
+ |
|
|
|
−1 |
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
+1 |
|
|
|
|
|
|
|
+1 |
|
|
|
|
|
|
|
|
||||||||||
|
|
|
|
|
|
|
( |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
− |
|
|
|
−1 |
|
|
|
|||||
|
|
∆ + ∆ |
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
∆ |
|
|
|
2 |
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
−1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
Ф |
− Ф |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
−1 |
) = 0 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
∆ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
−1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
2 |
|
|
|
|
|
|
−+1 + |
|
|
− |
|
|
Ф |
+1 |
− Ф |
|
|
|
− + |
|
−1 |
−−1 |
||||||||||||||||
|
|
|
|
( |
|
|
|
|
+1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
− |
|
|
|
|
|
|
|
||||||
|
∆ + ∆ |
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
∆ |
|
|
|
2 |
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
−1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Ф |
|
− Ф |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
−1 |
) = 0 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
∆ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
−1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
Граничные условия:
|
|
|
Ф |
|
|
|
= 1 |
|
|
|
|
||||
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|||
|
|
|
Ф |
|
= 1 |
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
Ф |
|
= − |
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
{ |
Ф |
|
|
= |
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
1 |
|
|
|
√ |
|
1 |
|
|
|
|
||
|
= ( |
|
+ |
|
|
) |
+ 1) |
||||||||
|
|
( |
|
||||||||||||
1 |
|
|
2 |
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
√ |
|
|
|
2 |
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
= ( |
|
+ |
|
( |
|
|
) + 1) + |
||||||||
|
|
|
|
||||||||||||
|
|
2 |
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
Текст программы и результаты метода Ньютона:
clear all close all clc warning off
constants; % инициализация необходимых физических констант
parameters; % инициализация шага на сетке, значения относительной погрешности к которой стремимся,
22
structure; % инициализация структуры, длины p и n областей, концентрации доноров и
акцепторов в данных областях,
%шага для подгонки, напряжения приложенного к диоду, расчет концентраций
%собственых носителей заряда, а также параметров необходимых для расчета
%подвижности электронов и дырок
grid_functions; % задание сетки координат, концентраций и подвижностей с разным шагом
для каждой области.
%загон конечной длины получившейся сетки в переменную I, а также
%инициализация матрицы разностей двух координат (x(2) - x(1)) и т.д.
normirovka; % нормировка координат, матрицы разностей, напряжения, концентрации,
подвижностей
prog_param; % значения для толщин линий и размера шрифта, название зашитое в переменной
% для будущих файлов
fi0(1)=log(N(1)/2+sqrt((N(1)/2)^2+1));% граничные условия для потенциала. В начале fi0(I)=-log(-N(I)/2+sqrt((N(I)/2)^2+1));% В конце
for i=2:I-1 % проход по всей сетке кроме граничных случаев, расчет электростатических
% потенциалов (далее электропотенциал) в ней fi0(i)=fi0(1)+(fi0(I)-fi0(1))/(x(I)-x(1))*(x(i)-x(1));
end
fi0=fi0'; % транспонирование
E=zeros(I,1); % инициализация матрицы единиц для напряженности электрического
% поля с I строк (кол-во точек в сетке) и 1 столбцом for j=1:length(U) % проход по всем заданным напряжениям
if j>1 % условие для задания начальных значений, если напряжение является не 1 по списку
for i=1:I
fi0(i)=fi0(i)+(fi0(i)-fi0(1))/(fi0(I)-fi0(1))*(U(j)-U(j-1)); % перерасчет электропотенциалов
end
Fn0=ones(I,1); % сброс прошлой матрицы для экспоненты квазиуровня ферми для электронов + задание начального значения
Fn0(I)=exp(-U(j)); % задание конечного значения для экспоненты квазиуровня ферми для электронов
Fp0=ones(I,1);
Fp0(I)=exp(U(j)); % повтор описанного выше, но уже для дырок
f0=[fi0; Fn0; Fp0]; % запись всех требующихся для нас параметров в одну
матрицу end
ct=0; % инициализация счетчика
23
err=100*Delta; % задание ошибки такой, чтоб начал работать цикл while err>Delta
ct=ct+1; % увеличение счетчика на единицу
if U(j)==0 % проверка для определения 0 ли напряжение по списку, и если 0 то
Fn0=ones(I,1); % инициализация матрицы для экспоненты квазиуровня ферми для электронов + задание начального значения
Fp0=ones(I,1); % тоже самое но для дырок
A=zeros(I,I); % инициализация матрицы Якоби
B=zeros(I,1); % инициализация матрицы значений потенциала
A(1,1)=1; % задание начального значения матрицы Якоби
B(1)=fi0(1)-log(N(1)/2+sqrt((N(1)/2)^2+1)); % задание начального значения матрицы значений потенциала
A(I,I)=1; % задание конечного значения матрицы Якоби
B(I)=fi0(I)+log(-N(I)/2+sqrt((N(I)/2)^2+1))-U(j); % задание конечного значения матрицы значений потенциала
for i=2:I-1 % проход по всем элементам матриц с расчетом каждой из них
A(i,i)=2/(dx(i)+dx(i-1))*(-1/dx(i)-1/dx(i-1))-Fn0(i)*exp(fi0(i))-
Fp0(i)*exp(-fi0(i));
A(i,i+1)=2/(dx(i)+dx(i-1))/dx(i); % расчет частной производной по
fi(i+1) - fi(i)
A(i,i-1)=2/(dx(i)+dx(i-1))/dx(i-1); % расчет частной производной по
fi(i) - fi(i-1)
B(i)=2/(dx(i)+dx(i-1))*((fi0(i+1)-fi0(i))/dx(i)-(fi0(i)-fi0(i-
1))/dx(i-1))-...
Fn0(i)*exp(fi0(i))+Fp0(i)*exp(-fi0(i))+N(i); % расчет частной производной по потенциалу для формулы градиента потенциала
end
fi=fi0-A^(-1)*B; % расчет разницы между начальными условиями и рассчитанными
% (-1 в степени А - обратная матрица)
clc % стирание из командного окна всей информации
Voltage=U(j)*V0 % денормирование напряжения, вывод на экран err=max(abs(fi-fi0))/max(abs(fi0)) % заменяем невязку на абсолютную
погрешность
ct % вывод счетчика итераций
fi0=fi0+(fi-fi0)./K(j); % расчет новых начальных условий else % если напряжение не 1 по списку - считаем по этой части
A=zeros(3*I,3*I); % инициализация матрицы Якоби сразу для
% электропотенциала и экспонент квазиуровня ферми
B=zeros(3*I,1); % инициализация матрицы значений функций (fi0, Fn0 и Fp0)
24
% задание начальных значений для
A(1,1)=1; % матрицы Якоби, отвечающей за fi0
B(1)=fi0(1)-log(N(1)/2+sqrt((N(1)/2)^2+1)); % матрицы значений для fi0 % задание конечных значений для:
A(I,I)=1; % матрицы Якоби, отвечающей за fi0
B(I)=fi0(I)+log(-N(I)/2+sqrt((N(I)/2)^2+1))-U(j); % матрицы значений для
fi0
% задание начальных значений для:
for i=2:I-1 % расчет значений для матрицы Якоби
%и расчет значений для матрицы значений функции
%частные производные для формулы градиента потенциала
A(i,i)=2/(dx(i)+dx(i-1))*(-1/dx(i)-1/dx(i-1))-Fn0(i)*exp(fi0(i))- Fp0(i)*exp(-fi0(i)); % расчет частной производной по fi(i)
A(i,i+1)=2/(dx(i)+dx(i-1))/dx(i); % расчет частной производной по
fi(i+1) - fi(i)
A(i,i-1)=2/(dx(i)+dx(i-1))/dx(i-1); % расчет частной производной по
fi(i) - fi(i-1)
A(i,I+i)=-exp(fi0(i)); % расчет частной производной по Fn0
A(i,2*I+i)=exp(-fi0(i)); % расчет частной производной по Fp0
B(i)=2/(dx(i)+dx(i-1))*((fi0(i+1)-fi0(i))/dx(i)-(fi0(i)-fi0(i-
1))/dx(i-1))-...
Fn0(i)*exp(fi0(i))+Fp0(i)*exp(-fi0(i))+N(i); % расчет значения
функции
%частные производные для формулы градиента экспоненты
%квазиуровня Ферми для электронов
end A(I+1,I+1)=1; B(I+1)=Fn0(1)-1; A(2*I,2*I)=1;
B(2*I)=Fn0(I)-exp(-U(j)); for i=2:I-1
A(I+i,I+i)=-(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2/dx(i)-...
(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fn(i)
A(I+i,I+i+1)=(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2/dx(i); %
расчет частной производной по Fn(i+1)-Fn(i) A(I+i,I+i-1)=(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2/dx(i-1); %
расчет частной производной по Fn(i)-Fn(i-1)
A(I+i,i)=mun(i)/2*exp(fi0(i))*(Fn0(i+1)-Fn0(i))/dx(i)-...
mun(i)/2*exp(fi0(i))*(Fn0(i)-Fn0(i-1))/dx(i-1); % расчет частной производной по fi(i)
25
A(I+i,i+1)=mun(i+1)/2*exp(fi0(i+1))*(Fn0(i+1)-Fn0(i))/dx(i); % расчет частной производной по fi(i+1)
A(I+i,i-1)=-mun(i-1)/2*exp(fi0(i-1))*(Fn0(i)-Fn0(i-1))/dx(i-1); %
расчет частной производной по fi(i-1) B(I+i)=(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2*(Fn0(i+1)-
Fn0(i))/dx(i)-...
(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2*(Fn0(i)-Fn0(i- 1))/dx(i-1); % расчет значения функции
%частные производные для формулы градиента экспоненты
%квазиуровня Ферми для электронов
end A(2*I+1,2*I+1)=1; B(2*I+1)=Fp0(1)-1; A(3*I,3*I)=1;
B(3*I)=Fp0(I)-exp(U(j)); for i=2:I-1
A(2*I+i,2*I+i)=-(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-
fi0(i)))/2/dx(i)-...
(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fp(i)
A(2*I+i,2*I+i+1)=(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(- fi0(i)))/2/dx(i); % расчет частной производной по Fp(i+1)-Fp(i)
A(2*I+i,2*I+i-1)=(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i- 1)))/2/dx(i-1); % расчет частной производной по Fp(i)-Fp(i-1)
A(2*I+i,i)=-mup(i)/2*exp(-fi0(i))*(Fp0(i+1)-Fp0(i))/dx(i)+...
mup(i)/2*exp(-fi0(i))*(Fp0(i)-Fp0(i-1))/dx(i-1); % расчет частной производной по fi(i)
A(2*I+i,i+1)=-mup(i+1)/2*exp(-fi0(i+1))*(Fp0(i+1)-Fp0(i))/dx(i); %
расчет частной производной по fi(i+1) A(2*I+i,i-1)=mup(i-1)/2*exp(-fi0(i-1))*(Fp0(i)-Fp0(i-1))/dx(i-1); %
расчет частной производной по fi(i-1) B(2*I+i)=(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-fi0(i)))/2*(Fp0(i+1)-
Fp0(i))/dx(i)-...
(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2*(Fp0(i)-Fp0(i- 1))/dx(i-1); % расчет значения функции
end
f=f0-A^(-1)*B; % расчет разницы между начальными условиями и рассчитанными
clc % очистка командного окна
Voltage=U(j)*V0 % денормирование напряжения, вывод на экран
26
err=max(abs(f-f0))/max(abs(f0)) % заменяем невязку на абсолютную
погрешность
ct % вывод счетчика итераций
f0=f0+(f-f0)./K(j); % расчет новых начальных условий fi0=f0(1:I);
Fn0=f0(I+1:2*I);
Fp0=f0(2*I+1:3*I); % замена предыдущих значений новыми
end
n=Fn0.*exp(fi0); % расчет концентрации электронов p=Fp0.*exp(-fi0); % расчет концентрации дырок for i=1:I-1
E(i)=-(fi0(i+1)-fi0(i))/dx(i); % расчет напряженности электрического поля
end
E(I)=E(I-1)+dx(I-1)/dx(I-2)*(E(I-1)-E(I-2)); % так как не можем посчитать последнее значение для
%поля ввиду отсутсвия следующего значения потенциала - принимаем
%его за последний рассчитанный
Ro=p-n+N;
end
for i=2:I-1 % так как плотность тока меняться внутри диода не будет -
% принимаем всю плотность тока за рассчитанную
Jn(i)=mun(i)*exp(fi0(i))*(Fn0(i+1)-Fn0(i-1))/(dx(i)+dx(i-1)); % расчет плотности тока электронов
Jp(i)=-mup(i)*exp(-fi0(i))*(Fp0(i+1)-Fp0(i-1))/(dx(i)+dx(i-1)); % расчет плотности ток дырок
end Jn(I)=Jn(I-1); Jp(I)=Jp(I-1); Jn(1)=Jn(2); Jp(1)=Jp(2);
J=Jn+Jp; % суммарная плотность тока
Current(j)=J(I); % запоминаем плотность тока для построения ВАХ graphics;% вывод графиков
end
figure % создание отдельного окна plot(U.*V0,-Current.*(J0*1e-4),'LineWidth',LW) % построение ВАХ xlabel('Voltage, V','FontSize',FS) % именуем ось абсцисс ylabel('Current density, A/cm^2','FontSize',FS) % именуем ось ординат grid on % включаем сетку на графике
print(gcf,'-djpeg','VAH') % сохраняем график ВАХ в виде файла джипег
Результаты расчетов:
27
Voltage = 0.8000
err = 2.8880e-08
Полученные графики:
Рис. 6. – Полученные программой графики.
Рис. 7. – Полученный программой график.
28
Текст программы и результаты метода Гуммеля:
clear all close all clc warning off
cconstants; % инициализация необходимых физических констант
parameters; % инициализация шага на сетке, значения относительной погрешности к которой стремимся,
structure; % инициализация структуры, длины p и n областей, концентрации доноров и акцепторов в данных областях,
%шага для подгонки, напряжения приложенного к диоду, расчет концентраций
%собственых носителей заряда, а также параметров необходимых для расчета
%подвижности электронов и дырок
grid_functions; % задание сетки координат, концентраций и подвижностей с разным шагом
для каждой области.
%загон конечной длины получившейся сетки в переменную I, а также
%инициализация матрицы разностей двух координат (x(2) - x(1)) и т.д.
normirovka; % нормировка координат, матрицы разностей, напряжения, концентрации,
подвижностей
prog_param; % значения для толщин линий и размера шрифта, название зашитое в переменной
% для будущих файлов
fi0(1)=log(N(1)/2+sqrt((N(1)/2)^2+1));% граничные условия для потенциала. В начале fi0(I)=-log(-N(I)/2+sqrt((N(I)/2)^2+1));% В конце
for i=2:I-1 %проход по всей сетке кроме граничных случаев, расчет электростатических
% потенциалов (далее электропотенциал) в ней fi0(i)=fi0(1)+(fi0(I)-fi0(1))/(x(I)-x(1))*(x(i)-x(1));
end
fi0=fi0';% транспонирование
E=zeros(I,1);% инициализация матрицы единиц для напряженности электрического
% поля с I строк (кол-во точек в сетке) и 1 столбцом for j=1:length(U) % проход по всем заданным напряжениям
if j>1 % условие для задания начальных значений, если напряжение является не 1 по списку
for i=1:I
fi0(i)=fi0(i)+(fi0(i)-fi0(1))/(fi0(I)-fi0(1))*(U(j)-U(j-1)); % перерасчет электропотенциалов
end
end
29
ct=0; % инициализация счетчика
err=100*Delta; % задание ошибки такой, чтоб начал работать цикл while err>Delta
ct=ct+1; % увеличение счетчика на единицу
if U(j)==0 % проверка для определения 0 ли напряжение по списку, и если 0 то
Fn=ones(I,1); % инициализация матрицы для экспоненты квазиуровня
% ферми для электронов + задание начального значения
Fp=ones(I,1); % тоже самое но для дырок else
A=zeros(I,I); % инициализация матрицы для частных производных Fn
B=zeros(I,1); % инициализация матрицы для значений функции
A(1,1)=1;
B(1)=1; % начальные условия
A(I,I)=1;
B(I)=exp(-U(j)); % конечные условия for i=2:I-1
A(i,i)=-(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2/dx(i)-...
(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2/dx(i-1); % расчет частной произврдная по Fn(i)
A(i,i+1)=(mun(i+1)*exp(fi0(i+1))+mun(i)*exp(fi0(i)))/2/dx(i); %
расчет частной производной по Fn(i+1)-Fn(i) A(i,i-1)=(mun(i)*exp(fi0(i))+mun(i-1)*exp(fi0(i-1)))/2/dx(i-1); %
расчет частной производной по Fn(i)-Fn(i-1) B(i)=0; % значение функции
end
Fn=A^(-1)*B; % решение уравнения
A=zeros(I,I);
B=zeros(I,1); % сброс матриц для использования в расчетах Fp
A(1,1)=1;
B(1)=1; % начальные условия
A(I,I)=1;
B(I)=exp(U(j)); % конечные условия for i=2:I-1
A(i,i)=-(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-fi0(i)))/2/dx(i)-...
(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2/dx(i-1); % расчет частной производной по Fp(i)
A(i,i+1)=(mup(i+1)*exp(-fi0(i+1))+mup(i)*exp(-fi0(i)))/2/dx(i); %
расчет частной производной по Fp(i+1)-Fp(i) A(i,i-1)=(mup(i)*exp(-fi0(i))+mup(i-1)*exp(-fi0(i-1)))/2/dx(i-1); %
расчет частной производной по Fp(i)-Fp(i-1)
30