[X-Shadow 0]

Здравствуйте. Столкнулся с проблемой при решении системы дифференциальных уравнений которая представлена на рисунке в приложении. Эта система содержит два интервала. Из за того что это система уравнений второго порядка и я не знаю значений для первых производных на границах (y'(-20) и y'(10)) то это задача с граничными условиями и в теории может быть решена с помощью функции bvp4c. Проблема в том как прописать эту систему уравнений вместе с условием непрерывности в районе точки 0 (снизу изображения). Я попытался сделать следующее (функция для передачи в bvp4c):

 

function [ dydx ] = tode( x , y )
dydx = zeros(2,1); % create zero array
dydx(1) = y(2);
if (x > 0)
    dydx(2) = b/a*((1+y(1))^(3/2)-1);
elseif (x <= 0)
    dydx(2) = 0;
end;
end

 

Однако эта функция не включает в себя условие с первыми производными. Подскажите пожалуйста как можно решить данную проблему? Спасибо.

[fider 0]

На вскидку:

Если эта система эквивалентна одному дифуру второго порядка, то двух дополнительных (граничных) условий должно хватить.

Похоже, что раз вторая производная на левом участке равна нулю, то там решение - прямая.

Непрерывность при (0) напоминает условие стыковки сплайнов.

А по сути - сегодня чуть позже попробую набросать решение (если получится).

Изменено 14 января, 2014 пользователем fider

[X-Shadow 0]

Вообще это задача посути есть моделирование электрического потенциала под затвором в МДП транзисторе. Участок где вторая производная равна нулю есть диэлектрик и там действительно согласно закону потенциал меняется линейно. В другой же части моделируется полупроводник. Из за аккумуляции носителей заряда возле границы полупроводник-диэлектрик (точка 0) потенциал там меняется не линейно. Условия y(-0)=y(+0) и a*y'(-0)=c*y'(+0) есть условия непрерывности потенциала и вектора электрической индукции соответственно. Большое спасибо, был бы очень благодарен за какой нибудь маленький примерчик

[X-Shadow 0]

Похоже что матлаб врятли сможет решить это встроенными функциями. Они все расчитаны на то что производная как и сама функция гладкие а следуя условию a*y'(-0)=c*y'(+0) выходит что первая производная кусочно гладкая. Поэтому врятли это можно решить встроенными функциями, скорее всего придется писать свой алгоритм

[fider 0]

Мquote name='X-Shadow' date='Jan 15 2014, 11:05' post='1226262']

Похоже что матлаб врятли сможет решить это встроенными функциями. Они все расчитаны на то что производная как и сама функция гладкие а следуя условию a*y'(-0)=c*y'(+0) выходит что первая производная кусочно гладкая. Поэтому врятли это можно решить встроенными функциями, скорее всего придется писать свой алгоритм

 

Мне тоже кажется, что лучше написать своими руками. С помощью встроенных функций все равно просто вряд ли получится.

Хотя конечно можно накрутить решение подгонкой (итерациями) на двух интервалах, чтобы в итоге они стыковались.

Или методом "стрельбы" или еще как-то.

 

И еще о задаче. Ведь речь идет о распределении потенциала по затвором по глубине канала?

Вроде-бы это рассматривалось. Можно, например, посмотреть google.1024.ru или еще где-нибудь.

 

Кстати, условие y(-0)=y(+0) в этой системе выполняется автоматически при y=0.

 

В то же время интересно. А это обязательно надо в МАТЛАБе? Это принципиально?

Изменено 15 января, 2014 пользователем fider

[X-Shadow 0]

Я задал этот вопрос на хелпе матлаба и там человек предложил решение http://www.mathworks.com.au/matlabcentral/...ential-equation. Опуская матлабовские функции этот мужик предлагает решить подобную задачу следующим образом (как я понял): разбить систему на 2 уравнения и решать их как будто 2 разных дифура с нач условиями. Значения первых производных не известны поэтому для начала приравнять обе к 0 (касательно правой части отрезка это верно т.к там она и равна нулю). Затем решить 2 уравнения отдельно друг от друга и проверить условие в районе точки 0. Если не сходятся то изменить начальное условие (увеличить или уменьшить начальное значение производной) и повторять пока не будут выполняться условия в районе точки 0. Все бы хорошо кроме одного нюанса - я точно знаю что производная в правом участке y'(10)=0 (ну или около нуля т.к в этом месте полупроводник имеет контакт с металлом и на границе т.к это оммический переход потенциалы должны быть равны а следовательно первая производная от функции потенциала тоже равна 0), это условие из физики. Т.е как бы мы это и не меняем. Но вот получается что второе уравнение никак не зависит от изменения параметров в первом. Фактически мы просто подгоняем решение первого уравнения чтобы оно сходилось с решением второго хотя на самом деле решение второго уравнения должно зависеть от граничного условия в точке х=-20 (в этой точке мы прикладываем потенциал U). Или может я чтото не так понял?

[fider 0]

Да.

Я примерно так и набросал, правда в МАТКАДе. (Прикрепляю). _____2.rtf_____2.rar

Только сначала я взял при х=-20 задал у=5 (в примере это z0), а потом условно-потолочно взял и задал руками y' (это z1) при х=-20.

Потом запустил как решение задачи Коши, то есть на интервале от -20 до 0 с этими начальными условиями, заданными при х=-20.

При этом использовал уравнение для диэлектрика.

В примере использовал функцию rkfixed (в МАТЛАБ должна быть аналогичная "ode45").

В результате получил распределение в диэлектрике и главное - при х=-0 получил значение y (это Rezult1(N1,1)) и y' (это Rezult1(N1,2)).

Здесь вообще Rezult - это массивы решения по точкам (здесь N1+1=100+1 точка решения).

 

После того как получил y и y' при х=-0, эти значения задал на границе для полупроводника как начальные.

То есть y(+0)=y(-0) и y'(+0)=(c/a)*y'(-0). И дальше продолжил решение для полупроводника от х=0 до х=10.

Подгонкой вручную (подгонял y' (это z1) при х=-20) получил в самой правой точке y=0.

 

Конечно, такой способ подгонки решения ненормальный. Но раз я неосторожно написал, что попробую - попробовал.

Все-таки это задача краевая и по-хорошему численно ее решать надо-бы наверное методом конечных разностей.

Но это надо еще повозиться. Но непонятно пока - зачем?

Ведь есть пакеты моделирования типа TCAD или хорошие примеры наработаны в flexPDE (прошлая ссылка).

 

По поводу условий сшивки решения на границе и на контакте. Примерно вроде бы так.

Электрическая индукция не должно претерпевать разрыв на границе двух сред.

А на контакте (в идеальном металле поле равно 0) касательная составляющая напряженности E конечно равна 0, с нормальной сложнее.

А потенциал, который связан через градиент с напряженностью поля может быть и не равен нулю?

 

[X-Shadow 0]

Большое спасибо за помощь!

 

Потенциал на границе х=10 должен быть равен нулю т.к там металл который подключен к земле. Производная тоже должна быть равна нулю т.к внутри металла электрического поля нет (рассматриваем идеальный случай).

 

Касательно почему не в пакетах TCAD. Я сначала хочу получить аналитический результат а потом сравнить его с симуляцией.

 

Сейчас пытаюсь реализовать это в MATLAB но пока как то туго. Собственно код:

 

global coef_A coef_B;

coef_A = 1.5477e+10;
coef_B = 9.7420e+03;

U = 1;              %значение y(-20e-7)
dy1 = 0;           %значение y'(-20e7)

err = 1e-1;       %допустимая ошибка значения y(10e-7)
end_val = 2;     %значение в точке y(10e-7)

xspan_neg = [-20e-7 : 1e-8 : -1e-8];      %интервалы
xspan_pos = [1e-8 : 1e-8 : 10e-7];

while end_val > err         %пока конечное значение больше ошибки
    
    y0_neg = [U dy1];                                                         %вектор потенциала и производной в точке -20у-7
    [X1,Y1] = ode45('neg_diff_eq',xspan_neg,y0_neg);          %вычисляем функцию y(x) и y'(x) на участке [-20e-7; 0]

    U1 = Y1(end,1);                                                             %значение y(+0), первое условие непрерывности
    dU1 = 0.5391*Y1(end,2);                                                %значение y'(+0), второй условие непрерывности

    y0_pos = [U1 dU1];                                                        %формируем вектор
    [X2,Y2] = ode45('pos_diff_eq',xspan_pos,y0_pos);           %вычисляем функцию y(x) и y'(x) на участке [0; 10e-7]

    end_val = abs(Y2(end,1));                                              %полученное значение значение y(10e-7)
    dy1 = dy1 - 1e4;                                                            %уменьшаем значение первой производной
    
end


X = [X1; X2]; Y = [Y1; Y2];                                                %склеиваем графики

subplot(2,1,1);                                                                  %рисуем
plot(X,Y(:,1));
subplot(2,1,2);
plot(X,Y(:,2));

 

А сами функции выглядят так:

 

function [ dydx ] = neg_diff_eq( x , y )

dydx = zeros(2,1); % create zero array
dydx(1) = y(2);
dydx(2) = 0;
            
end

function [ dydx ] = pos_diff_eq( x , y )

global coef_A coef_B;

dydx = zeros(2,1); % create zero array
dydx(1) = y(2);

dydx(2) = coef_A*(1+coef_B*y(1))^(3/2)-1;
            
end

 

Но в итоге решение почемуто не сходится. Матлаб просто зависает пытаясь это решить

Изменено 17 января, 2014 пользователем X-Shadow