Интеграл и график в любой версии Pascal, Maple. Что-либо другое

[Dimafania]

Интеграл и график в любой версии Pascal, Maple. Что-либо другое

Здравствуйте. Буду премного благодарен тем, кто решит поставленную передо мной задачу. Мне необходимо решить определенный интеграл вида A(x,y)B(x)dx и построить его график. Проблема заключается в том, что функция B(x) у меня задана в виде набора значений x. Справиться с подобным интегралом пока что не в состоянии. Заранее спасибо

 

[Vladimir_S]

Справиться с подобным интегралом пока что не в состоянии. Заранее спасибо

Вам просто посочувствовать, или как? Потому что для помощи в решении задачи нам бы увидеть её. Прям как есть целиком. С функциями, таблицами, полной формулировкой задания и всем таким прочим.
Ну а пока - сочувствуем. Душой с Вами.

 

[Dimafania]

Вам просто посочувствовать, или как? Потому что для помощи в решении задачи нам бы увидеть её. Прям как есть целиком. С функциями, таблицами, полной формулировкой задания и всем таким прочим.
Ну а пока - сочувствуем. Душой с Вами.

Здравствуйте, Владимир.
Спасибо за ваш ответ.
К своему сообщению прикладываю два файла. В первом приведена формула с интегралом, во втором приведены значения функции g(w).
Могу привести вам значения различных констант из формулы (как например постоянная Больцмана Кв), но поскольку меня больше интересует алгоритм, а вам данные значения интереса не представляют (вы можете их заменить единицей например), приводить их, с вашего разрешения, не буду.
Суть задачи заключается в следующем. В прикрепленном файле .txt находятся значения плотности фононных состояний g(w), или, иными словами, количество фононов обладающих определенной энергией. Используя данную плотность состояния можно определить теплоемкость материала. Формула связи приведена на рисунке. Граничными значениями w являются соответственно граничные значения файла (w0 = 0.730, wmax = 8.487). Последовательно подставляя различные значения температуры Т можно получить график теплоемкости, либо просто набор значений. Если бы функция g(w) задавалась бы какой-нибудь функцией, то проблем бы не возникло. Однако в моей случае это набор значений. Понятно, что интеграл можно заменить суммой произведения всех значений функции на шаг дискретизации, однако в моей случае существует проблема понимания того, как написать код в котором часть интеграла задается функцией, а вторая часть - суммой.

 

[Dimafania]

Приложения к сообщению

 

[Vladimir_S]

Здравствуйте, Владимир. Спасибо за ваш ответ.

Да не за что.

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

Вообще-то я сам физик, поэтому с постоянными Планка и Больцмана проблем нет. А вот значение величины V привести бы стоило.

Понятно, что интеграл можно заменить суммой произведения всех значений функции на шаг дискретизации, однако в моей случае существует проблема понимания того, как написать код в котором часть интеграла задается функцией, а вторая часть - суммой.

Ну... да, правильно - и при этом я напрочь не вижу проблем. Вот, расписал алгоритм:
Посмотреть вложение 226767
Программу сами напишете, или как?

 

[Dimafania]

Большое вам спасибо, Владимир.
Программированием, по правде говоря, начал заниматься буквально пару дней назад, но программу постараюсь написать самостоятельно (в дальнейшем это, не сомневаюсь, пригодится). В случае если где-то опять встряну в тупик, то напишу вам.
За объем я буду брать единицу.
В приведенном вами алгоритме все понятно, лишь один момент попрошу прокомментировать:
Почему в итоговом интеграле в конце дельта w вы берете вычитанием (w[i+1]-w) без деления на 2, а в функциях F и g идет суммирование с делением на 2? Судя по моей логике везде должна быть такая форма (w[i+1]-w)

 

[Vladimir_S]

Большое вам спасибо, Владимир.
Программированием, по правде говоря, начал заниматься буквально пару дней назад, но программу постараюсь написать самостоятельно (в дальнейшем это, не сомневаюсь, пригодится). В случае если где-то опять встряну в тупик, то напишу вам.
За объем я буду брать единицу.
В приведенном вами алгоритме все понятно, лишь один момент попрошу прокомментировать:
Почему в итоговом интеграле в конце дельта w вы берете вычитанием (w[i+1]-w) без деления на 2, а в функциях F и g идет суммирование с делением на 2? Судя по моей логике везде должна быть такая форма (w[i+1]-w)

Дык - геометрия требует. Площадь трапеции, знаете ли. Вот:

 

[Dimafania]

Спасибо еще раз, Владимир. Ваша помощь бесценна)
Насчет аппроксимации я думал использовать не трапецию,а обычные прямоугольники (хотя понятно, что точность при этом будет ниже).

 

[Dimafania]

Владимир, вынужден снова обратиться к вам за помощью.
К сожалению написанная мною программа выдает какую-то ерунду.
Ниже приведу вам текст, может быть вы увидите какую-нибудь глупость.
Перед тем как приведу вам текст, хочу обратить внимание, что файл я переименовал в Zr.txt, а также оставил там только значения g(w), удалив значения w. Значения хотел подсчитать от 1 до 10К (для начала). Текст программы:
program Specific_Heat;
const
dE : Comp = 0.033E12;
h : Comp = 6.582119E-16;
K : Real =8.6173324E-5;
var
t : text;
x:array [1..239] of real;
i, Tem, j :integer;
F,s,Ei : Comp;
begin
assign(t, 'Zr.txt');
reset(t);
i:=0;
while not eof (t) do
begin
i:=i+1;
readln(t,x);
writeln(x);
end;
s:=0;
for i:=1 to 239 do
begin
Ei:=0.730E12+(i*dE-dE/2);
F:=Ei*Ei/sqr((exp(h*Ei/2*K*Tem)-exp(-h*Ei/2*K*Tem))/2);
s:=s+F*x;
end;
s:=(dE*s*h*h)/(4*K*Tem*Tem);
j:=1;
for j:=1 to 10 do
begin
Tem:=j;
writeln(s);
end;
end.

 

[Dimafania]

Comp изменил на Extended - не помогло...

 

[Vladimir_S]

Владимир, вынужден снова обратиться к вам за помощью.

Пожалуйста, чем, как говорится, могу.

К сожалению написанная мною программа выдает какую-то ерунду.

Это нормально. Так всегда бывает у начинающих программистов.

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

Сначала я напишу свой вариант, а Вы посмотрите.

Перед тем как приведу вам текст, хочу обратить внимание, что файл я переименовал в Zr.txt, а также оставил там только значения g(w), удалив значения w.

Вот уж совсем не могу понять - зачем? Ведь шкала частот - неравномерная!
P.S. Какие-то странные у Вас значения K и ħ. Это в каких единицах?

 

[Dimafania]

Шкалу частот я принял равномерно увеличивающейся с шагом 0.033E12.
Постоянную Больцмана брал в эВ*К^-1
Планка эВ*с
Подбирал размерности таким образом, чтобы гиперболический синус был безразмерным.

 

[Vladimir_S]

Шкалу частот я принял равномерно увеличивающейся с шагом 0.033E12.
Постоянную Больцмана брал в эВ*К^-1
Планка эВ*с
Подбирал размерности таким образом, чтобы гиперболический синус был безразмерным.

Значит, так, сударь.
Написал я программку и прокрутил.
Вот:

Const
 h=6.582119E-16;
 K=8.6173324E-5;
 N=239;

Var
 t:text;
 w,g:Array [1..239] of real;
 i,Tem:integer;

Function Sh(r:real):real;
begin
 Sh:=(Exp(r)-Exp(-r))/2;
end;

Function F(v,T:real):real;
begin
 F:=K*Sqr(h*v/(K*T*2))/Sqr(Sh(h*v/(K*T*2)));
end;

Function Cv(T:real):Real;
var
 Sum:Real;
 j:Integer;
begin
 Sum:=0;
 for j:=1 to N-1 do
  Sum:=Sum+F((w[j+1]+w[j])/2,T)*(g[j+1]+g[j])/2*(w[j+1]-w[j]);
 Cv:=Sum;
end;

begin
 Assign(t,'D:\g.txt');
 Reset(t);
 Readln(t);
 for i:=1 to N do readln(t,w[i],g[i]);
 Close(t);
 Assign(t,'D:\res.txt');
 Rewrite(t);
 for Tem:=1 to 10 do
  Writeln(t,Tem:3,'     ',Cv(1.0*Tem));
 Close(t)
End.

Пояснения: исходный файл g.txt я оставил, как был. Выходные данные отправил в файл Res.txt. Оба файла у меня расположены на диске D: (эту уж мне так удобно; Вы поставьте по-своему). В результате получил практически константу, что немудрено: аргумент гиперболического синуса имеет порядок 10^-10, а это значит, что гиперболический синус можно заменить самим аргументом, в результате чего в подынтегральном выражении функция F становится просто равной К, чего быть не должно.
Вывод: что-то напортачено с единицами в константах. Сейчас попробую разобраться.

 

[Vladimir_S]

А, нет, прошу прощения - отбой! Константы поставлены правильно, проверил. Это просто я чудак на букву "м" - забыл домножить частоту на 10^12. Если это сделать, получается вполне разумно. Вот исправленный текст:

Const
 h=6.582119E-16;
 K=8.6173324E-5;
 N=239;

Var
 t:text;
 w,g:Array [1..239] of real;
 i,Tem:integer;

Function Sh(r:real):real;
begin
 Sh:=(Exp(r)-Exp(-r))/2;
end;

Function F(v,T:real):real;
begin
 F:=K*Sqr(h*v/(K*T*2))/Sqr(Sh(h*v/(K*T*2)));
end;

Function Cv(T:real):Real;
var
 Sum:Real;
 j:Integer;
begin
 Sum:=0;
 for j:=1 to N-1 do
  Sum:=Sum+F((w[j+1]+w[j])/2,T)*(g[j+1]+g[j])/2*(w[j+1]-w[j]);
 Cv:=Sum;
end;

begin
 Assign(t,'D:\g.txt');
 Reset(t);
 Readln(t);
 for i:=1 to N do
  begin
   readln(t,w[i],g[i]);
   w[i]:=w[i]*1E12;
  end;
 Close(t);
 Assign(t,'D:\res.txt');
 Rewrite(t);
 for Tem:=1 to 10 do
  Writeln(t,Tem:3,'     ',Cv(1.0*Tem));
 Close(t)
End.

 

[Dimafania]

Спасибо, Владимир. Сейчас буду списывать вашу программу и проверять) Интересно, что я попробовал на пальцах проверить свою программу (только интеграл взял полегче, чтоб вручную подсчитать), и вроде бы все работает. А с изначальным интегралом брехня получается=)

 

[Dimafania]

Sum:=Sum+F((w[j+1]+w[j])/2,T)*(g[j+1]+g[j])/2*(w[j+1]-w[j]);

Владимир, а в этой строчке множитель (w[j+1]-w[j]) не уходит случайно к двойке в знаменатель? Скобками не надо выделять?

 

[Vladimir_S]

Владимир, а в этой строчке множитель (w[j+1]-w[j]) не уходит случайно к двойке в знаменатель? Скобками не надо выделять?

Нет. В Паскале (да и, по-моему, в любом ЯП) умножение и деление имеют равный приоритет, и цепочка вида a*b/c/d*e означает: a умножается на b, результат делится на c, результат делится на d, результат умножается на e.

 

[Vladimir_S]

Для пущей наглядности: программа

Begin
 Writeln(2*3/2*3);
 Readln
End.

выдает в результате 9, а не 1. Можете проверить.

 

[Dimafania]

Еще раз спасибо, Владимир. Программа ваша работает прекрасно. Буду разбираться с полученными цифрами, а точнее степенями.

 

[Dimafania]

Writeln(t,Tem:3,' ',Cv(1.0*Tem));

Добрый день, Владимир. Не могли бы вы пояснить эту строчку? Точнее член Cv(1.0*Tem)? А еще точнее вот этот момент 1.0*Tem? Значение теплоемкости емкости умножается на значение температуры?