Свертка трехмерной плотности распределения вероятности.

@dobryasha · Регистрация: 04.02.2013

Студворк — интернет-сервис помощи студентам

Здравствуйте, уважаемые форумчане.

Мне необходимо решить следующую задачу:

Мне известна трехмерная плотность распределения вероятностей (ПРВ):

$https://www.cyberforum.ru/cgi-bin/latex.cgi?W({u}_{1},{u}_{2},{u}_{3})=1,353*\frac{{u}_{1}*{u}_{2}*{u}_{3}}{{\sigma}^{6}}*{exp}^{-\frac{1,394*\left({{u}_{1}}^{2}+{{u}_{3}}^{2} \right)+1,826*{{u}_{2}}^{2}}{2{\sigma}^{2}}}*{I}_{0}\left(-0,845*\frac{{u}_{1}*{u}_{2}}{{\sigma}^{2}}\right)*{I}_{0}\left(0,334*\frac{{u}_{1}*{u}_{3}}{{\sigma}^{2}}\right)*{I}_{0}\left(-0,845*\frac{{u}_{2}*{u}_{3}}{{\sigma}^{2}}\right)$

Мне известны также выражения для $https://www.cyberforum.ru/cgi-bin/latex.cgi?{u}_{1},{u}_{2},{u}_{3}$ . Они отличаются друг от друга временными моментами. Т.е. выражение для $https://www.cyberforum.ru/cgi-bin/latex.cgi?{u}$ записывается для трех моментов времени $https://www.cyberforum.ru/cgi-bin/latex.cgi?{t}_{1},{t}_{2},{t}_{3}$ . И мы исследуем зависимость этих трех отсчетов.

Так вот мне нужно перейти от трехмерной ПРВ к одномерной. Это можно сделать с помощью двух операций свертки. Например, с помощью первой свертки прийти к выражению вида $https://www.cyberforum.ru/cgi-bin/latex.cgi?W({u}_{1},{u}_{2})$ , а с помощью второй к $https://www.cyberforum.ru/cgi-bin/latex.cgi?W({u}_{1})$ . Я знаю, что в Матлабе есть функция $https://www.cyberforum.ru/cgi-bin/latex.cgi?y = conv(x, h)$ , но она легко используется для числовых последовательностей.

А вот как с помощью Матлаба свернуть нужное мне выражение я пока не смекаю. Подскажите пожалуйста, если у Вас есть мысли на эту тему.
Заранее спасибо.

@dobryasha · 28.02.2013, 14:36 **[ТС]**

А тут не получится такого $https://www.cyberforum.ru/cgi-bin/latex.cgi?\sqrt{Q1(i)-\sqrt{U1}}$ ?

Зосима · 28.02.2013, 16:10

нет

потому, как я сначала нахожу корень разности $https://www.cyberforum.ru/cgi-bin/latex.cgi?\sqrt{O1(i)-U1}$
А потом уже заменяю U1 на корень!

@dobryasha · 28.02.2013, 18:02 **[ТС]**

Ну и хорошо) Еще раз огромное спасибо!!

Добавлено через 1 час 12 минут
Етишкин пистолет!!
Последний график не строится...Ввожу разные интервалы, пустота

Зосима · 28.02.2013, 19:23

Подправил предел, ведь U1 должно меняться от 0 до О1:

Matlab M
1
U1 = linspace(0, O1(i), 20);

Но там все-равно в Z получаются все NaN-ы, которые plot игнорирует...

Добавлено через 18 минут
Попробуй это:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
clear, clc
O1 = linspace(1e-15, 1e-9, 20);
Z = zeros(size(O1));
for i = 1:length(O1)
    U1 = linspace(1e-15, O1(i), 20);
    Q1 = sqrt(O1(i)-U1);
    U1 = sqrt(U1);
    z = funW(U1)./(2*U1).*funW(Q1)./(2*Q1);
    z(isnan(z)) = 0; % обнуляем NaNы
    Z(i) = trapz( U1, z );
end
plot(O1,Z)

Поясню позже

Пояснение: в последний момент, когда нужно было быстро собираться и бежать на заводскую маршрутку меня осенило

1. В якобиане u1 в знаменателе, поэтому при u1=0 естессно полечаем юринду, поэтому я заменил ноль на близкое к нему значение 1e-15

2. Подинтегральную ф-цию рассчитал в отдельную переменную z, чтобы при интегрировании убрать NaN-ы, вернее обнулить: z(isnan(z))=0; и потом уже считать интеграл trapz.
Впосделствии можно даже убрать точкозапятую, чтобы результаты отображались на экране.

3. Методом научного тыка было установлено, что при верхнем пределе O1 большим примерно 1e-5 все значения z получаются NaN-ами, поэтому взял 1e-9, но возможно надо будет с ним поиграться, чтоб поиметь красивый график

Но как ты догадалась, запустить и посмотреть результат я не успел

Поэтому с огромным нетерпением жду от тебя ответа, шо ж там таки выходит!

@dobryasha · 01.03.2013, 07:58 **[ТС]**

Я вчера после работы пришла домой и рухнула спать)) Так что пока не запустила новый кодец. Сейчас запущу

Добавлено через 11 минут
Запустила. Опять пустой экран, шоб его

С 1e-9 точно надо поиграться. Сейчас поиграюсь)

Добавлено через 10 минут
Я ввела 1e-4. График получился такой: от 1e-5 до 1e-4 ноль. От 0 до 1e-5 наклонная прямая, уходящая в бесконечность (до 2*10^18)

Добавлено через 36 секунд
График выложить не могу, т.к. матлаб на ноуте, а на нем нет интернета.

Добавлено через 34 секунды
А на рабочем компе usb разьем отключен

Добавлено через 3 минуты
А если посмотреть еще область меньше нуля, то получается треугольничек с вершиной в точке примерно [5e-6;1.8*10^18]

Зосима · 01.03.2013, 08:41

спинным мозгом чую: якобиан всю малину портит!

Слух, а попробуй убрать корни, якобианы и оставить чисто свертку:

U1 = linspase...
z = funW(U1).*funW(O1(i)-U1)
z(...) = 0;
trapz ... и т.д.

@dobryasha · 01.03.2013, 15:54 **[ТС]**

эх, друг мой...ездила к преподу, опять перестановка...возвращаюсь к первой программе (без сумматоров и квадраторов)
функция теперь выглядит так:

Matlab M
1
2
3
4
5
6
7
8
9
10
tau1 = 1.268e-4;
tau2 = 1.268e-4;
tau3 = 1.268e-4;
fm = 0.738*(U1(i)*(U2(j)-U1(i)).*(U3-U2(j)))/(tau1*tau2*tau3*sigma^6).* ...
            exp( -( 0.123*(U1(i)/tau1).^2 + 0.123*((U3-U2(j))/tau3).^2 + 0.349*((U2(j)-U1(i))/tau2).^2 )/(2*sigma^2) );
        
        % вторая часть подинтегральной ф-ции
        b1 = -0.192*U1(i)*(U2(j)-U1(i))/(tau1*tau2*sigma^2);
        b2 = 0.095*U1(i)*(U3-U2(j))/(tau1*tau3*sigma^2);
        b3 = -0.192*(U2(j)-U1(i)).*(U3-U2(j))/(tau2*tau3*sigma^2);

При этом график ведет себя как ошалелый

А всего-то, поменяли значения констант!

@dobryasha · 02.03.2013, 13:31 **[ТС]**

График получается вот такой:

@dobryasha · 02.03.2013, 13:35 **[ТС]**

Я поняла, что дело в экспоненте...Как только хотя бы один из коэффициентов заменяешь, например 1.394 на 0.123, график улетает вниз.

@dobryasha · 02.03.2013, 13:36 **[ТС]**

А начало-то у него хорошее:

@dobryasha · 03.03.2013, 10:40 **[ТС]**

Нет, не хотя бы один коэффициент в степени экспоненты, а именно первый, перед u1/tau1 все портит, гад.

@dobryasha · 03.03.2013, 10:54 **[ТС]**

Вот, при U1 = linspace(0, 1.55e-3, 100) график выходит такой:

@dobryasha · 03.03.2013, 17:23 **[ТС]**

Т.е в определенный момент график начинает ускоренно стремитсья в минус бесконечность. А если коэффициент при u1/tau1 в степени экспоненты сделать например не 0.123, а 0.523, то график строится нормальным...Блин, как так-то!!

Добавлено через 1 час 26 минут
Да даже при 0.223 график нормальный строится

Добавлено через 5 часов 0 минут
В общем так. Ввела я импульсню характеристику фильтра. Графики получились нормальные.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
function W = funRELEY(U1)
 
sigma = 1;
tau1 = 0.028537;
tau2 = 0.102558;
tau3 = 0.220953;
Dtau = 1/(tau1*tau2*tau3);
 
W = zeros(size(U1));
for i = 1:length(U1)
    U2 = linspace(0,U1(i),70);
    F = zeros(size(U2));
    for j = 1:length(U2)
        U3 = linspace(0,U1(i)-U2(j),70);
        fm = 0.738*(U1(i)*(U2(j)-U1(i)).*(U3-U2(j)))/(tau1*tau2*tau3*sigma^6).* ...
            exp( -( 0.123*(U1(i)/tau1).^2 + 0.123*((U3-U2(j))/tau3).^2 + 0.349*((U2(j)-U1(i))/tau2).^2 )/(2*sigma^2) );
        
        % вторая часть подинтегральной ф-ции
        b1 = -0.192*U1(i)*(U2(j)-U1(i))/(tau1*tau2*sigma^2);
        b2 = 0.095*U1(i)*(U3-U2(j))/(tau1*tau3*sigma^2);
        b3 = -0.192*(U2(j)-U1(i)).*(U3-U2(j))/(tau2*tau3*sigma^2);
        S = 0;
        for k = 0:5
            % накапливаем сумму
            S = S + (-1).^k.*(1+(k>0)).*besseli(k,b1).*besseli(k,b2).*besseli(k,b3);
            % (1+(k>0)) - это символ Кронекера: =1 при k=0 и =2 при k>0
        end
        FdU3 = Dtau.*fm.*S.*(U2(j)>U3);
        % считаем интеграл по dU3
        F(j) = trapz(U3, FdU3);
    end
    % считаем bнтеграл по dU2
    W(i) = trapz(U2, F);
end

@dobryasha · 03.03.2013, 17:24 **[ТС]**

Вот для одномерной ПРВ:

Matlab M
1
2
3
4
5
6
7
8
U1 = linspace(0.5e-15, 0.7, 70);
W = funRELEY(U1);
hPlot = plot(U1, W, 'r');
set( hPlot, 'LineWidth', 3 );
title('W(u_1)','fontsize',20)
xlabel('u_1','fontsize',20)
ylabel('W','fontsize',20)
grid on

@dobryasha · 03.03.2013, 17:26 **[ТС]**

Это после квадратора:

Matlab M
1
2
3
4
5
6
7
U1 = linspace(0.5e-15, 0.3, 100);
K = funRELEY(sqrt(U1))./(2*sqrt(U1));
hPlot = plot(U1, K, 'r');
set( hPlot, 'LineWidth', 3 );
xlabel('k_1','fontsize',20)
ylabel('W(k_1)','fontsize',20)
grid on

@dobryasha · 03.03.2013, 17:27 **[ТС]**

После сумматора:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
clear, clc
O1 = linspace(1e-15, 0.7, 70);
Z = zeros(size(O1));
for i = 1:length(O1)
    U1 = linspace(1e-15, O1(i), 70);
    U1 = sqrt(U1);
    Q1 = sqrt(O1(i)-U1);
    z = funRELEY(U1)./(2*U1).*funRELEY(Q1)./(2*Q1);
    z(isnan(z)) = 0; % обнуляем NaNы
    Z(i) = trapz( U1, z );
end
hPlot = plot(O1, Z, 'r');
set( hPlot, 'LineWidth', 3 );
xlabel('w_1','fontsize',20)
ylabel('W(w_1)','fontsize',20)
grid on

@dobryasha · 04.03.2013, 11:12 **[ТС]**

А теперь еще осталось построить график ПРВ после блока извлечения квадратного корня.
Т.е. теперь нужно произвести замену x1 = sqrt(u1) на y1 = x1^2 плюс якобиан перехода 2y1. Пока не пойму, как из последнего кода перейти к нужному выражению

Думаю

Т.е. у нас было f(U1), потом f(sqrt(U1)), потом интеграл от f(Q1-sqrt(U1)) по dU1 , а теперь нужно f(Q1^2)*2Q1 (якобиан)

Добавлено через 17 часов 9 минут
Зосима, помоги пожалуйста, как время будет

@dobryasha · 04.03.2013, 21:26 **[ТС]**

У меня новая беда. Теперь препод хочет, чтобы осталась зависимость не от u1, а от u3. Я изменяю код, график слетает с катушек. Что я делаю не так?

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
function W = funnRELEY(U3)
 
sigma = 1;
tau1 = 0.028537;
tau2 = 0.102558;
tau3 = 0.220953;
Dtau = 1/(tau1*tau2*tau3);
 
W = zeros(size(U3));
for i = 1:length(U3)
    U1 = linspace(0,U3(i),20);
    F = zeros(size(U1));
    for j = 1:length(U1)
        U2 = linspace(0,U3(i)-U1(j),20);
        fm = 0.738*(U1(j)*(U2-U1(j)).*(U3(i)-U2))/(tau1*tau2*tau3*sigma^6).* ...
            exp( -( 0.123*(U1(j)/tau1).^2 + 0.123*((U3(i)-U2)/tau3).^2 + 0.349*((U2-U1(j))/tau2).^2 )/(2*sigma^2) );
        
        % вторая часть подинтегральной ф-ции
        b1 = -0.192*U1(j)*(U2-U1(j))/(tau1*tau2*sigma^2);
        b2 = 0.095*U1(j)*(U3(i)-U2)/(tau1*tau3*sigma^2);
        b3 = -0.192*(U2-U1(j)).*(U3(i)-U2)/(tau2*tau3*sigma^2);
        S = 0;
        for k = 0:5
            % накапливаем сумму
            S = S + (-1).^k.*(1+(k>0)).*besseli(k,b1).*besseli(k,b2).*besseli(k,b3);
            % (1+(k>0)) - это символ Кронекера: =1 при k=0 и =2 при k>0
        end
        FdU2 = Dtau.*fm.*S.*(U1(j)>U2);
        % считаем интеграл по dU3
        F(j) = trapz(U2, FdU2);
    end
    % считаем bнтеграл по dU2
    W(i) = trapz(U1, F);
end

Matlab M
1
2
3
4
5
6
7
8
9
10
U3 = linspace(0.5e-15, 7, 70);
W = funnRELEY(U3);
h = U3(2)-U3(1);
L = sum(W*h);
hPlot = plot(U3, W, 'r');
set( hPlot, 'LineWidth', 3 );
title('W(u_3)','fontsize',20)
xlabel('u_3','fontsize',20)
ylabel('W','fontsize',20)
grid on

@dobryasha · 04.03.2013, 21:29 **[ТС]**

Зосима, помоги пожааааалуйста, как время будет

Зосима · 05.03.2013, 10:46

Прости, лапка

с тырнетом были перебои и на работе чуток завал

Но у меня для тебя есть хорошая новость!

Свертка трехмерной плотности распределения вероятности.

Причесал программу, пришлось очень внимательно заменить везде в теле функции U1(i) на U3(i), U2(j) на U1(j) и U3 на U1 (только благодаря индексам i,j и старой ф-ции разобрался

) Держи код функции:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
function W = funnRELEY(U3)
 
sigma = 1;
tau1 = 0.028537;
tau2 = 0.102558;
tau3 = 0.220953;
Dtau = 1/(tau1*tau2*tau3);
 
W = zeros(size(U3));
for i = 1:length(U3)
    U1 = linspace(0,U3(i),20);
    F = zeros(size(U1));
    for j = 1:length(U1)
        U2 = linspace(0,U3(i)-U1(j),20);
        fm = 0.738*(U3(i)*(U1(j)-U3(i)).*(U2-U1(j)))/(tau1*tau2*tau3*sigma^6).* ...
            exp( -( 0.123*(U3(i)/tau1).^2 + 0.123*((U2-U1(j))/tau3).^2 + 0.349*((U1(j)-U3(i))/tau2).^2 )/(2*sigma^2) );
        
        % вторая часть подинтегральной ф-ции
        b1 = -0.192*U3(i)*(U1(j)-U3(i))/(tau1*tau2*sigma^2);
        b2 = 0.095*U3(i)*(U2-U1(j))/(tau1*tau3*sigma^2);
        b3 = -0.192*(U1(j)-U3(i)).*(U2-U1(j))/(tau2*tau3*sigma^2);
        S = 0;
        for k = 0:5
            % накапливаем сумму
            S = S + (-1).^k.*(1+(k>0)).*besseli(k,b1).*besseli(k,b2).*besseli(k,b3);
            % (1+(k>0)) - это символ Кронекера: =1 при k=0 и =2 при k>0
        end
        FdU2 = Dtau.*fm.*S.*(U1(j)>U2);
        % считаем интеграл по dU2
        F(j) = trapz(U2, FdU2);
    end
    % считаем bнтеграл по dU1
    W(i) = trapz(U1, F);
end

И программка расчета:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
clear, clc
U3 = linspace(0.5e-15, 1, 70);
W = funnRELEY(U3);
W(isnan(W))=0;
L = trapz(U3, W)
hPlot = plot(U3, W, 'r');
set( hPlot, 'LineWidth', 3 );
title('W(u_3)','fontsize',20)
xlabel('u_3','fontsize',20)
ylabel('W','fontsize',20)
grid on

Тут я только уменьшил верхний предел с 7 до 1, чтоб наглядней было

ну и чтоб интеграл считался, обнулил NaN-ы

По поводу обнуления NaN-ов - это не подгонка! Мы их просто исключаем из интегральной суммы, чтоб они ее не портили!

Новые блоги и статьи Все статьи Все блоги /
Саморегулирующийся социальный контракт для сервера cross-section. Hrethgir 14.08.2026 С кодом конечно таких глубоких размышлений пока не было, впрочем я уже привык к алгоритмизации. Суть предмета записи: снова в диалоге с нейросетью (я взял пока себе ник для учётки админа - Rector). . . .	Часы электронные Uhbif79 12.08.2026 Выкладываю программу часов. Программа позволяет: 1. Использовать системное время и дату, 2. Есть возможность вводить время и дату вручную. 3. Реализованы 2 будильника: начало и конец рабочего дня. . . .	Часы с будильником на основе класса QLCDNumber Uhbif79 12.08.2026 Всем добрый день, выкладываю программу часов с будильником на основе класса QLCDNumber. Здесь я пробовал самостоятельно создавал классы, впервые столкнулся с видимостью переменной одного класса из. . .	Установка MinGW GCC 16.2 и CMake 8Observer8 10.08.2026 VK Видео: https:/ / vkvideo. ru/ video-240781534_456239017 YouTube: eY5-5PyI9NM Текстовая версия
Неделя из жизни имитационной модели склада: мои кривые руки растут, откуда надо anaschu 10.08.2026 Неделя из жизни имитационной модели склада: как я почти написал неправильную логику и что с этим делать Работаю сейчас над учебно-рабочим проектом: строю в AnyLogic имитационную модель процессов. . .	Калькулятор для расчета родства russiannick 07.08.2026 1. Задача: Создать калькулятор для расчета родства. Родственных связей существует 8 ступеней, такие как: p - отец P - мать q - муж Q - жена b - брат B - сестра s - сын S - дочь	Мир по моей воле kumehtar 07.08.2026 Когда-то кажется, что всё просто. Ты весь такой светлый. Причиняешь добро. Борешься за справедливость в этом тёмном мире. Потом начинаешь замечать одну неприятную вещь. Почти каждый хороший. . .	Кредитный калькулятор Maks 05.08.2026 Решение задачи по прикладной информатике средствами 1С. Задача: Напишите приложение-калькулятор, которое помогает рассчитывать параметры кредита для аннуитетного и дифференцированного видов. . .

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101

	Свертка трехмерной плотности распределения вероятности. 04.02.2013, 11:54. Показов 15879. Ответов 163 Метки нет (Все метки) Здравствуйте, уважаемые форумчане. Мне необходимо решить следующую задачу: Мне известна трехмерная плотность распределения вероятностей (ПРВ): $https://www.cyberforum.ru/cgi-bin/latex.cgi?W({u}_{1},{u}_{2},{u}_{3})=1,353\frac{{u}_{1}{u}_{2}{u}_{3}}{{\sigma}^{6}}{exp}^{-\frac{1,394\left({{u}_{1}}^{2}+{{u}_{3}}^{2} \right)+1,826{{u}_{2}}^{2}}{2{\sigma}^{2}}}{I}_{0}\left(-0,845\frac{{u}_{1}{u}_{2}}{{\sigma}^{2}}\right){I}_{0}\left(0,334\frac{{u}_{1}{u}_{3}}{{\sigma}^{2}}\right){I}_{0}\left(-0,845\frac{{u}_{2}*{u}_{3}}{{\sigma}^{2}}\right)$ Мне известны также выражения для $https://www.cyberforum.ru/cgi-bin/latex.cgi?{u}_{1},{u}_{2},{u}_{3}$ . Они отличаются друг от друга временными моментами. Т.е. выражение для $https://www.cyberforum.ru/cgi-bin/latex.cgi?{u}$ записывается для трех моментов времени $https://www.cyberforum.ru/cgi-bin/latex.cgi?{t}_{1},{t}_{2},{t}_{3}$ . И мы исследуем зависимость этих трех отсчетов. Так вот мне нужно перейти от трехмерной ПРВ к одномерной. Это можно сделать с помощью двух операций свертки. Например, с помощью первой свертки прийти к выражению вида $https://www.cyberforum.ru/cgi-bin/latex.cgi?W({u}_{1},{u}_{2})$ , а с помощью второй к $https://www.cyberforum.ru/cgi-bin/latex.cgi?W({u}_{1})$ . Я знаю, что в Матлабе есть функция $https://www.cyberforum.ru/cgi-bin/latex.cgi?y = conv(x, h)$ , но она легко используется для числовых последовательностей. А вот как с помощью Матлаба свернуть нужное мне выражение я пока не смекаю. Подскажите пожалуйста, если у Вас есть мысли на эту тему. Заранее спасибо. 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	28.02.2013, 14:36 [ТС]
	А тут не получится такого $https://www.cyberforum.ru/cgi-bin/latex.cgi?\sqrt{Q1(i)-\sqrt{U1}}$ ? 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	28.02.2013, 16:10
	нет потому, как я сначала нахожу корень разности $https://www.cyberforum.ru/cgi-bin/latex.cgi?\sqrt{O1(i)-U1}$ А потом уже заменяю U1 на корень! 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	28.02.2013, 18:02 [ТС]
	Ну и хорошо) Еще раз огромное спасибо!! Добавлено через 1 час 12 минут Етишкин пистолет!! Последний график не строится...Ввожу разные интервалы, пустота 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	01.03.2013, 07:58 [ТС]
	Я вчера после работы пришла домой и рухнула спать)) Так что пока не запустила новый кодец. Сейчас запущу Добавлено через 11 минут Запустила. Опять пустой экран, шоб его С 1e-9 точно надо поиграться. Сейчас поиграюсь) Добавлено через 10 минут Я ввела 1e-4. График получился такой: от 1e-5 до 1e-4 ноль. От 0 до 1e-5 наклонная прямая, уходящая в бесконечность (до 210^18) Добавлено через 36 секунд* График выложить не могу, т.к. матлаб на ноуте, а на нем нет интернета. Добавлено через 34 секунды А на рабочем компе usb разьем отключен Добавлено через 3 минуты А если посмотреть еще область меньше нуля, то получается треугольничек с вершиной в точке примерно [5e-6;1.8*10^18] 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	01.03.2013, 08:41
	спинным мозгом чую: якобиан всю малину портит! Слух, а попробуй убрать корни, якобианы и оставить чисто свертку: U1 = linspase... z = funW(U1).*funW(O1(i)-U1) z(...) = 0; trapz ... и т.д. 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	02.03.2013, 13:31 [ТС]
	График получается вот такой: Миниатюры 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	02.03.2013, 13:35 [ТС]
	Я поняла, что дело в экспоненте...Как только хотя бы один из коэффициентов заменяешь, например 1.394 на 0.123, график улетает вниз. 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	02.03.2013, 13:36 [ТС]
	А начало-то у него хорошее: Миниатюры 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	03.03.2013, 10:40 [ТС]
	Нет, не хотя бы один коэффициент в степени экспоненты, а именно первый, перед u1/tau1 все портит, гад. 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	03.03.2013, 10:54 [ТС]
	Вот, при U1 = linspace(0, 1.55e-3, 100) график выходит такой: Миниатюры 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	04.03.2013, 11:12 [ТС]
	А теперь еще осталось построить график ПРВ после блока извлечения квадратного корня. Т.е. теперь нужно произвести замену x1 = sqrt(u1) на y1 = x1^2 плюс якобиан перехода 2y1. Пока не пойму, как из последнего кода перейти к нужному выражению Думаю Т.е. у нас было f(U1), потом f(sqrt(U1)), потом интеграл от f(Q1-sqrt(U1)) по dU1 , а теперь нужно f(Q1^2)2Q1 (якобиан) Добавлено через 17 часов 9 минут* Зосима, помоги пожалуйста, как время будет 0

@dobryasha 0 / 0 / 0 Регистрация: 04.02.2013 Сообщений: 101
	04.03.2013, 21:29 [ТС]
	Зосима, помоги пожааааалуйста, как время будет 0