Интерполяция кусочнозаданной функции

@SBoris · Регистрация: 26.07.2012

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

Есть задача:
дан набор значений С1 = f(t) таблица значений.

Matlab M

if x<=2
  C1 = C0+Cexp*ka/ke*(1-exp(-ke*t));
else
  C1 = C0+Cexp*ka/ke*(exp(-ke*(t-tc))-exp(-ke*t));
end

Параметры для оптимизации ka,ke,tc.

Как это протянуть в cftool?
У меня идеи кончились. Отдельно - идёт, а вот вместе для определения tc...
Сейчас думаю попробовать сосорудить функцию с аргументами типа f(ka,ke,tc) и подсунуть в cftool. Правда, я не уверен сработает ли.
Если бы возможно было соорудить конструкцию оператора ветвления в одну строку кода.... Но, не знаю как.

Спасибо за ответы и соображения.

Зосима · 13.11.2012, 19:36

Попробуй создать файл-функцию fun(ka,ke,tc), и запихнуть ее в fittype - кажется это единственная возможность

Добавлено через 49 минут
Кстать! Ф-цию можно записать приспособив для этого ф-цию хевисайда! =)
C = C1 + (C2 -C1)*heaviside(x-2);
Где вместо C1 и C2 нужно записать выражения для x<2 и x>2 соответственно.
Единственное, пишу с телефона и могу ошибиться в имени ф-ции хевисайда(непростая у него фамилия), глянь help или гугл спроси.

Смысл заключается в следующем: когда x<2 ф-ция хевисайда равна 0 и C = C1, а когда x>2 ф-ция хевисайда равна 1 и C = C1 + C2 - C1 = C2, т.е. то что надо!

Добавлено через 19 секунд
Кстать! Ф-цию можно записать приспособив для этого ф-цию хевисайда! =)
C = C1 + (C2 -C1)*heaviside(x-2);
Где вместо C1 и C2 нужно записать выражения для x<2 и x>2 соответственно.
Единственное, пишу с телефона и могу ошибиться в имени ф-ции хевисайда(непростая у него фамилия), глянь help или гугл спроси.

Смысл заключается в следующем: когда x<2 ф-ция хевисайда равна 0 и C = C1, а когда x>2 ф-ция хевисайда равна 1 и C = C1 + C2 - C1 = C2, т.е. то что надо!

@SBoris · 14.11.2012, 20:30 **[ТС]**

Да, thnks,

- это есть гуд. НО, возникла проблема (или я заработался

)

Matlab M

function to_fit = C(t,Co,Cexp,ka,ke,tc)
 to_fit = Co+Cexp*ka/ke*(1-exp(-ke*t)) + (Co+Cexp*ka/ke*(exp(-ke*(t-tc))-exp(-ke*t))-Co+Cexp*ka/ke*(1-exp(-ke*t)))*heaviside(t-tc);
end
 
Co = 10.6; Cexp = 376; ka = 0.13; ke = 0.15; tc = 2;
t = 0.1:0.1:15;
ff = to_fit(t,Co,Cexp,ka,ke,tc);
 
filenameTKmean = sprintf('/home/boris/Documents/TK_data.csv');
dataTKmean = csvread (filenameTKmean);
tTKmean = dataTKmean(:,1); C_Ni_TKmean = dataTKmean(:,2);
 
plot(t,ff,'linewidth',3,tTKmean, C_Ni_TKmean, 'o');

дает неожиданный график (рис.1 в атаче)

А правильная отрисовка получается токо:

Matlab M

function to_fit = C(t,Co,Cexp,ka,ke,tc)
 to_fit = Co+Cexp*ka/ke*(1-exp(-ke*t)) + (Co+Cexp*ka/ke*(exp(-ke*(t-tc))-exp(-ke*t))-Co+Cexp*ka/ke*(1-exp(-ke*t)))*heaviside(t-tc);
end
 
Co = 10.6; Cexp = 376; ka = 0.13; ke = 0.15; tc = 2;
ff=to_fit(0,Co,Cexp,ka,ke,tc);
tt = 0;
for t = 0.1:0.1:15;
ff = [ff; to_fit(t,Co,Cexp,ka,ke,tc)];
tt = [tt; t];
end;
 
filenameTKmean = sprintf('/home/boris/Documents/TK_data.csv');
dataTKmean = csvread (filenameTKmean);
tTKmean = dataTKmean(:,1); C_Ni_TKmean = dataTKmean(:,2);
 
plot(tt,ff,'linewidth',2,tTKmean, C_Ni_TKmean, 'o');

ИМХО, этот подарок ещё лучше условного оператора в фиттинге. Где я заглючил - хоть убейте не вижу.

Зосима · 14.11.2012, 20:43

А можешь файлик TK_data.csv выложить в zip архиве, чтоб я мог на досуге поиграться.
*пока что я в коде не разобрался, как-то мудрено вышло (явно не хватает 50мл вдохновения ;D )

@SBoris · 14.11.2012, 21:02 **[ТС]**

В cftool получил

Мда... что-то не ладится с работой.

@SBoris · 14.11.2012, 21:10 **[ТС]**

Сообщение от Зосима

А можешь файлик TK_data.csv выложить в zip архиве, чтоб я мог на досуге поиграться.
*пока что я в коде не разобрался, как-то мудрено вышло (явно не хватает 50мл вдохновения

)

И я о том же. Были бы ближе, дринькнули бы.

Зосима · 15.11.2012, 11:50

Кстать, я не очень понял задание

Если ты делаешь интерполяцию, зачем тебе рассчитывать исходную функцию?
Ну да ладно

я ее оставил для красоты и сравнения.
Вот что вышло:

Matlab M

clear
clc
 
% Записываем исходную функцию
Cfun = @(t,tc,C0,Cexp,ka,ke) C0+Cexp*ka/ke*(1-exp(-ke*t)) +...
    heaviside(t-tc).*( (C0+Cexp*ka/ke*(exp(-ke*(t-tc))-exp(-ke*t)))...
    -(C0+Cexp*ka/ke*(1-exp(-ke*t))) );
 
% загружаем данные
TK = load('TK_data.csv');
xdata = TK(:,1);
ydata = TK(:,2);
 
% вводим параметры функции
C0 = 10.6; Cexp = 376; ka = 0.13; ke = 0.15; tc = 2;
t = 0:0.01:25;
 
% рассчитываем исходную функцию
y = Cfun(t,tc,C0,Cexp,ka,ke);
 
% создаем интерполяционную модель на базе данной функции
myfit = fittype('C0+Cexp*ka/ke*(1-exp(-ke*x)) + heaviside(x-tc).*( (C0+Cexp*ka/ke*(exp(-ke*(x-tc))-exp(-ke*x)))-(C0+Cexp*ka/ke*(1-exp(-ke*x))) )   ')
% производим интерполяцию данных
Ffit = fit(xdata, ydata, myfit)
% рассчитываем значения интерполированной ф-ции
yf = Ffit(t);
 
% Рисуем все
plot(xdata,ydata,'or',t,y,':g',t,yf,'-b')
legend('Точки данных', 'Исходная ф-ция', 'Интерополированная ф-ция')
 
 
 
% Проверочка. 
% Считаем отклонения от данных
 
y1 = Cfun(xdata,tc-0.01,C0,Cexp,ka,ke); % отсчеты исходной ф-ции
y2 = Ffit(xdata); % отсчеты интерполированной ф-ции
 
r1 = sum((ydata - y1).^2)/1000 % сумма квадратов разностей исходной ф-ции и данных
r2 = sum((ydata - y2).^2)/1000 % сумма квадратов разностей интерп-ой ф-ции и данных
if r1>r2
    disp('Интерполированная ф-ция ближе к данным')
elseif r1<r2
    disp('Исходная ф-ция ближе к данным')
else
    disp('Ф-ции равноудалены')
end

У матлабовской интерполяции есть особенность:
При каждом запуске результат разный и иногда совсем далекий от реальности,
поэтому не пугайся, а просто еще раз запусти программу!

(Иногда приходится и перезапускать матлаб!)

Будут вопросы - задавай!

@SBoris · 15.11.2012, 15:15 **[ТС]**

Если ты делаешь интерполяцию, зачем тебе рассчитывать исходную функцию?

Для самопроверки правильности отображения функции. См.выше - была проблема с отрисовкой функции с условным оператором.

Это для чего в 37ой строке? Оно кардинально влияет на результат.

Matlab M

tc-0.01

Зосима · 15.11.2012, 15:27

В массиве даных есть точка x=2. И как раз в этой точке происходит скачек ф-ции хевисайда, там получается 2 значения, неопределенность, матлаб не может сообразить, что делать и выдает результат NaN (нечисло).
При построении графика матлаб отбрасывет эту точку, поэтому там все более-менее гладко, а при точном расчете ф-ции получается ошибка, поэтому я чуток схитрил: сдвинул скачек на 0.01, что не должно сильно влиять на результат. (можно написать еще меньше - будет точнее).

Новые блоги и статьи
В чем отличие между INNER JOIN и OUTER JOIN bytestream 22.01.2025 В современных базах данных информация часто распределена между множеством взаимосвязанных таблиц, что делает операции объединения JOIN неотъемлемой частью работы с SQL. Эти операции позволяют. . .	Как сделать первую букву заглавной в JavaScript bytestream 22.01.2025 JavaScript предоставляет разработчикам множество инструментов для эффективной работы с текстовыми данными. Одной из часто встречающихся задач при обработке строк является преобразование первой буквы. . .	Что такое Big O нотация и алгоритмическая сложность bytestream 22.01.2025 Введение в алгоритмическую сложность В мире разработки программного обеспечения эффективность алгоритмов играет crucial роль в создании качественных приложений. Алгоритмическая сложность. . .	Как решать конфликты слияния (merge) в Git bytestream 22.01.2025 Конфликты слияния в системе контроля версий Git возникают в ситуациях, когда две или более ветки разработки содержат несовместимые изменения в одних и тех же участках кода. Эти конфликты представляют. . .	Как использовать регулярные выражения bytestream 22.01.2025 Регулярные выражения представляют собой мощный инструмент для работы с текстовыми данными, который позволяет осуществлять поиск, проверку и манипуляцию строками на основе определенных шаблонов. Этот. . .	Как выйти из Vim bytestream 22.01.2025 Vim (Vi IMproved) представляет собой один из самых влиятельных текстовых редакторов в истории компьютерной индустрии, эволюционировавший из своего предшественника Vi, созданного Биллом Джоем в 1976. . .
NoSQL базы данных: что это такое и какие существуют bytestream 22.01.2025 В современную эпоху цифровой трансформации объемы данных растут экспоненциально, создавая новые вызовы для традиционных систем управления базами данных. NoSQL (Not Only SQL) представляет собой. . .	Обновление исследования от команды MCM (январь 2025 г.) Programma_Boinc 22.01.2025 Обновление исследования от команды MCM (январь 2025 г. ) Мы продолжаем изучать молекулярные сигнатуры, связанные с раком легких, с текущим фокусом на GCM1, факторе транскрипции, участвующем в. . .	Как работать с Kafka в Go (Golang) bytestream 22.01.2025 Apache Kafka представляет собой распределенную платформу потоковой передачи данных, которая произвела революцию в области обработки событий и интеграции микросервисов. Эта система, изначально. . .	Как использовать RabbitMQ в Go (Golang) bytestream 22.01.2025 RabbitMQ представляет собой надежный и широко используемый брокер сообщений, который играет ключевую роль в построении современных распределенных систем и микросервисной архитектуры. В основе работы. . .	Как преобразовать список списков в простой список в Python bytestream 22.01.2025 При работе с Python разработчики часто сталкиваются с необходимостью обработки сложных структур данных, среди которых особое место занимают вложенные списки. Эти структуры представляют собой списки,. . .	Что такое GUID / UUID и как их создать bytestream 22.01.2025 В мире разработки программного обеспечения существует постоянная потребность в уникальной идентификации объектов, записей и ресурсов. Эта задача становится особенно актуальной в распределенных. . .

Зосима 5243 / 3571 / 379 Регистрация: 02.04.2012 Сообщений: 6,474 Записей в блоге: 17
	13.11.2012, 19:36	2
	Попробуй создать файл-функцию fun(ka,ke,tc), и запихнуть ее в fittype - кажется это единственная возможность Добавлено через 49 минут Кстать! Ф-цию можно записать приспособив для этого ф-цию хевисайда! =) C = C1 + (C2 -C1)heaviside(x-2); Где вместо C1 и C2 нужно записать выражения для x<2 и x>2 соответственно. Единственное, пишу с телефона и могу ошибиться в имени ф-ции хевисайда(непростая у него фамилия), глянь help или гугл спроси. Смысл заключается в следующем: когда x<2 ф-ция хевисайда равна 0 и C = C1, а когда x>2 ф-ция хевисайда равна 1 и C = C1 + C2 - C1 = C2, т.е. то что надо! Добавлено через 19 секунд* Кстать! Ф-цию можно записать приспособив для этого ф-цию хевисайда! =) C = C1 + (C2 -C1)*heaviside(x-2); Где вместо C1 и C2 нужно записать выражения для x<2 и x>2 соответственно. Единственное, пишу с телефона и могу ошибиться в имени ф-ции хевисайда(непростая у него фамилия), глянь help или гугл спроси. Смысл заключается в следующем: когда x<2 ф-ция хевисайда равна 0 и C = C1, а когда x>2 ф-ция хевисайда равна 1 и C = C1 + C2 - C1 = C2, т.е. то что надо! 1

Зосима 5243 / 3571 / 379 Регистрация: 02.04.2012 Сообщений: 6,474 Записей в блоге: 17
	14.11.2012, 20:43	4
	А можешь файлик TK_data.csv выложить в zip архиве, чтоб я мог на досуге поиграться. *пока что я в коде не разобрался, как-то мудрено вышло (явно не хватает 50мл вдохновения ;D ) 1

@SBoris 0 / 0 / 0 Регистрация: 26.07.2012 Сообщений: 20
	14.11.2012, 21:02 [ТС]	5
	В cftool получил Мда... что-то не ладится с работой. 0

Зосима 5243 / 3571 / 379 Регистрация: 02.04.2012 Сообщений: 6,474 Записей в блоге: 17
	15.11.2012, 15:27	9
	В массиве даных есть точка x=2. И как раз в этой точке происходит скачек ф-ции хевисайда, там получается 2 значения, неопределенность, матлаб не может сообразить, что делать и выдает результат NaN (нечисло). При построении графика матлаб отбрасывет эту точку, поэтому там все более-менее гладко, а при точном расчете ф-ции получается ошибка, поэтому я чуток схитрил: сдвинул скачек на 0.01, что не должно сильно влиять на результат. (можно написать еще меньше - будет точнее). 1