Решение системы уравнений (маткад)

[GAA]

alcoholist, спасибо! Я без особых оснований предположил, что $\Delta U_k>0$. Если $\Delta U_k < 0$, то интеграл будет сходиться.
И графики рисуются.

-= Добавлено через 2 часа =-

Грубо (по 9 точкам)

Вложение:

UkJk.PNG [ 5.13 Кб | Просмотров: 2556 ]

Вполне возможно, где-то мог допустить ошибку.

g:=0.15: p := 66; A := 3; B := 34; e0:= 8.85e-12; bi:= 5e-6; p := 66;
dk :=4^(1/3)*bi^(1/3)*U^(2/3)*e0^(1/3)*(1+g)^(1/3)/jk^(1/3);
II := A*p*Int(exp(B/(2*U/p/dk*(x-dk)/dk)), x=0..dk) - ln(1+1/g);

Точки --- нули II.

[optimden]

К сожалению, речь идёт о физике плазмы - $\delta U_k$ - падение напряжения в прикатодной области, оно там только положительным может быть.
Хотя... Сейчас снова глянул в пособие. Там кажется замена неправильно сделана. Изначально в степени экспоненты стоит:

$$\exp\left(\frac{-B}{(\frac{E0}p)((\frac{x}{d_k})-1)}\right)$$

$$E0=\frac{-2\delta U_k}{d_k}$$

То есть по-хорошему степень экспоненты благодаря минусу, который появится при подстановке, станет положительной. И тогда задача действительно будет разрешима.
Большое спасибо!
P.S. А в какой программе вы построили график?

[optimden]

В общем вопрос о нерешаемости данной задачи отпал. Были найдены опечатки в формулах. Спасибо всем, кто помогал.
Теперь возник другой вопрос, по части построения графика решения в матлабе:

Используется синтаксис Matlab M

i=1;
for jk=0.00001:0.01:10
    dk=@(deltaUk)((4)^(1/3)*bi^(1/3)*deltaUk^(2/3)*eps0^(1/3)*(y+1)^(1/3))/(jk^(1/3));
    f2=@(deltaUk) (A_const.*p.*int(exp((B_const.*p.*(dk.^2))./(2.*deltaUk.*(x-dk))),x,0,dk)-log(1+(1./y)));
    Uk=fzero(f2,[0 3000]);
    Mass1(1,i)=jk;
    Mass2(2,i)=Uk;
    i=i+1;
end
disp(Mass1);
 

На это матлаб выдаёт ошибку:

FZERO cannot continue because user supplied function_handle ==>
@(deltaUk)(A_const.*p.*int(exp((B_const.*p.*(dk.^2))./(2.*deltaUk.*(x-dk))),x,0,dk)-log(1+(1./y)))
failed with the error below.

Undefined function or method 'power' for input arguments of type 'function_handle'.

Error in ==> Collision at 91
    Uk=fzero(f2,[0 3000]);

Подскажите, пожалуйста, в чём может быть причина? Вроде и $d_k$, и $f2$ определены, как функции $\delta U_k$.

[GAA]

Эта ошибка связана с оплошностью при записи вызова функции: вместо вызова функции dk(deltaUk) Вы возводите в степень handle (dk).

optimden в сообщении #1387063 писал(а):

f2=@(deltaUk) (A_const.*p.*int(exp((B_const.*p.*(dk.^2))./(2.*deltaUk.*(x-dk))),x,0,dk)-log(1+(1./y)));

Не только в степень возводите handle. Исправить нужно все вызова функции. После исправления этой оплошности появятся другие сообщения об ошибках.

[optimden]

Я правильно понимаю, что надо просто заменить $d_k$ на $d_k(deltaUk)$ в теле функции $f2$? Просто такая замена не устранила ошибку. Насколько понимаю, при использовании в функции $f2$ просто переменной $d_k$ , матлаб воспринимает её не как функцию от $deltaUk$, а как указание на эту функцию.

[GAA]

Да, всюду в теле f2 надо заменить dk на dk(deltaUk). После этого сообщение об ошибке должно измениться.

[optimden]

Собственно сообщение об ошибке изменилось на такое:

??? Error using ==> fzero at 245
FZERO cannot continue because user supplied function_handle ==>
@(deltaUk)(A_const.*p.*int(exp((B_const.*p.*(dk(deltaUk).^2))./(2.*deltaUk.*(x-(dk(deltaUk))))),x,0,dk(deltaUk))-log(1+(1./y)))
failed with the error below.

Undefined function or variable 'x'.

Error in ==> Collision at 91
    Uk=fzero(f2,[0.00001 3000]);

Суть ошибки в том, что функция fzero не может работать с функциями, определёнными через анонимные?

[GAA]

Грубо говоря, в matlab определённый интеграл можно вычислить численно или символьно (Symbolic Math Toolbox).

Для численного вычисления в старых версиях использовались, например, функции quad (на основе формулы Симпсона) quadl (adaptive Lobatto quadrature), а в новых версиях — integral.

Для символьного интегрирования используется функция int. Подынтегральное выражение (первый входной аргумент int) должно быть символьным. Если Вы планируете использовать символьное интегрирование, то перед циклом можно создать символьную переменную x. Например, при помощи
syms x
или при помощи конструктора символьного объекта
x =sym('x')
В таком случае выражение, которое задаёт подынтегральную функцию (первый входной аргумент функции int), будет создано как символьное выражение.

Результат, возвращаемый функцией int — символьное выражение, которое, возможно, содержит имена констант или функций. Для преобразования такого выражения в число в представлении арифметики переменной точности (Variable-Precision Arithmetic) используется функция vpa. (При помощи vpa не будет ещё получено число с плавающей точкой двойной точности.) [Вместо int + vpa можно попробовать воспользоваться vpaintegral.]

Для преобразования в Double можно воспользоваться функцией double (конструктором объекта типа double array).

(После данных исправлений, если все константы заданы, возможно заработает вызов функции f2. Перед использованием этой функции в функции fzero я бы проверил работоспособность f2.)

[optimden]

Ругается как на

Используется синтаксис Matlab M

syms x

, так и на

Используется синтаксис Matlab M

x=syms('x')

. Я программу пишу не в comman window, а в отдельном m-файле. Может в этом дело?

[GAA]

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

Возможно не установлен Symbolic Math Toolbox. Попробуйте так проверить. Если следующее вернёт 1, то установлен, иначе --- нет.
>> v = ver; any(strcmp('Symbolic Math Toolbox', {v.Name}))

[optimden]

Да, действительно не установлен. Отличный повод обновится до последней версии матлаба.

[GAA]

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

$$J = pAd_k \int_0^1 \exp \left( \frac {Bpd_k}{2\Delta U_k (u-1)} \right) du.$$

Для численного интегрирования не нужен Symbolic Math Toolbox. Например, можно воспользоваться функцией integral. Первым входным аргументом этой функции является handle подынтегрального выражения. Обозначим переменную, хранящую этот handle, через f. Проверим, будет ли работать функция для векторного аргумента deltaUk при фиксированном значении jk.
Текст для matlab R2013b ниже — идея реализации. (Возможно, не самая разумная идея реализации.) Значения констант и выражения надо проверить.

Используется синтаксис Matlab M

y =0.145; p = 66; A = 3; B = 34; eps0 = 8.85e-12; bi = 5.18e-6;

jk = 0.1;

dk=@(deltaUk)(4)^(1/3)*bi^(1/3)*deltaUk.^(2/3)*eps0.^(1/3)*(y+1).^(1/3)/(jk^(1/3));
f = @(u, deltaUk) exp(B*p*dk(deltaUk)./(2* deltaUk*(u-1)));
J = @(deltaUk) integral( @(u) f(u, deltaUk), 0,1, 'ArrayValued',true);
f2 = @(deltaUk) p*A*dk(deltaUk).*J(deltaUk) - log(1+1/y);
deltaUk = 10000:100:100000;
plot(deltaUk, f2(deltaUk));
xlabel('\Delta U_k'); ylabel('f2');
Uk=fzero(f2,[1 10000000])
 

Получающийся рисунок:

Вложение:

Комментарий к файлу: F2(Delta Uk) matlab 2013b

F2(deltaUk).png [ 4.83 Кб | Просмотров: 2482 ]

Редактирование: убрал лишние скобки в исходнике.

[optimden]

Спасибо большое, всё работает!