Scilab для кухонных физиков

[CASTRO]

" но чем так  плох fortran2003?"

Может и не плох. Но приемственности уже нет. Никто не будет на нём пепеписывать cernlib, geant3... Да и переносимость у фортрана плохая. ЦЕРН перешел на 64 битный SLC5 уж год как, а cernlib так под 64 бита закомпилить и не удалось. Хотя, про отдельных шаманов ходят легенды...

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

Факт заключается в том, что физика частиц, которая начиналась на фортране в 70х, сейчас с фортрана сходит. Ключевые пакеты переписаны на C++ (GEANT4 сменил GEANT3, PYTHIA 8 скоро окончательно заменит PYTHIA6, ROOT победил PAW...). Если эксперименты конца 90х-начала 2000х еще обрабатывались на фортране, то сейчас уже все на С++.

[ssp123]

По большому счету, языков программирования очень много. Как и всяких разных задач.  Под каждую задачу можно найти язык, оптимальный для ее решения. Под некоторые задачи даже разрабатывается отдельный язык.  С/С++ заточен на создание операционных систем. Паскаль хорош для обучения школьников. Бэйсик- язык для бейсичных маньяков. Хаскель- для фунционального программирования- решения логических задач, искусственного интеллекта и т.п.
Матлабоподобные языки хороши тем, что в них реализована матричная арифметика, они просты для освоения, и интерпретаторы - можно пользоваться как каклькулятором, и не нужно при этом ждать, пока программа скомпилируется. Кроме того,  они позволяют транслировать свой код в фортран или с, если уж кому-то этого захочется.

Даны ускорения g на различной высоте h
Из этой таблицы ускорений нужно составить уравнение  для скорости
и решить задачу  - какова скорость у поверхности, то есть на уровне h=0 ?

Это Вы можете продемонстрировать своим Scilab ?

Да легко. Только тема создана и для Вас в том числе, я когда увидел Вашу битву с калькулятором в одной руке, с карандашом - в другой, я понял- что такую тему создать просто необходимо.

Поставьте ROOT. Это фактически интерпретатор C++ c  обширными библиотеками классов. Хотите -матрицы перемножайте, хотите - 4-векторы. Гистограмы, фиты (в самом общем случае), красивые картинки, деревья, фурье-анализ, геометрия... Это стандартный пакет обработки данных, используемый в физике высоких энергий. Качать тут

http://root.cern.ch/drupal/

А самое главное, если библиотеки ROOT можно использовать при компилляции программ (во всяком случае, с gcc)

Ага, спасибо. Только вряд ли пригодится, я ушел из науки лет 20 назад, а в условиях кухни- достаточно  скайлаба. Или purebasic'а - меня прикалывает он своей реактивной скоростью, сравнимой с сишной и ассемблерной.

[ssp123]

Вектора задаются так:
-->a=[1, 3, 5, 7]
Это вектор- строка, элементы разделены запятыми. В качестве элементов могут быть числа, логические значения (T или F, истина или ложь) и даже строки.
Вектор-столбец задается так:
-->b=[2; 4; 6; 8]

Матрицы  задаются таким образом;
m=[1,2,4,8;  1,3,9,27;  1,4,16,64;  1,5,25,125]

то есть, матрица есть вектор столбец, состоящий из векторов - строк
если есть набор векторов-строк v1,v2,v3 то матрицу можно задать так
m=[v1;v2;v3]

[ssp123]

Виктор Вашкевич, вот Ваша задача.

G = 6.672d-11
M = 5.976d+24
R = 6.376370d+6

v=0 // Начальная скорость, м/с
dt=1 // Отрезок времени, через который производим вычисления, с
r=400000000 // Расстояние до центра земли, м
cnt=0 // счетчик итераций

while (r>R)
  g=G*M/(r*r)
  r0=r
  r=r-(v*dt + g*dt*dt*0.5)
  v=v+g*dt
  cnt=cnt+1
end

disp(cnt,"Итераций выполнено=")
disp(v, " Скорость на поверхности= ")
disp(g, " Ускорение на поверхности= ")
disp(r0, " Расстояние до центра Земли= ")


Вставляйте ее в свой скайлаб, и мучайте.

У меня получилось все как по классической науке:
Итераций выполнено=   
 
    444621. 
 
  Скорость на поверхности=   
 
    11093.754 
 
  Ускорение на поверхности=   
 
    9.7866945 
 
  Расстояние до центра Земли=   
 
    6382859.7 
 

[ssp123]

А чой-то я ерундой-то занялся? Есть литература по скайлабу по русски, есть.
http://teacher.dn-ua.com/Math/Scilab/Scilab.html
http://www.csa.ru/~zebra/my_scilab/

[БЫК]

Еще вот
http://linux.armd.ru/common/img/uploaded/files/ScilabBook.pdf

[Вашкевич Виктор]

Виктор Вашкевич, вот Ваша задача.

G = 6.672d-11
M = 5.976d+24
R = 6.376370d+6

v=0 // Начальная скорость, м/с
dt=1 // Отрезок времени, через который производим вычисления, с
r=400000000 // Расстояние до центра земли, м
cnt=0 // счетчик итераций

while (r>R)
  g=G*M/(r*r)
  r0=r
  r=r-(v*dt + g*dt*dt*0.5)
  v=v+g*dt
  cnt=cnt+1
end

disp(cnt,"Итераций выполнено=")
disp(v, " Скорость на поверхности= ")
disp(g, " Ускорение на поверхности= ")
disp(r0, " Расстояние до центра Земли= ")


Вставляйте ее в свой скайлаб, и мучайте.

У меня получилось все как по классической науке:
Итераций выполнено=   
 
    444621. 
 
  Скорость на поверхности=   
 
    11093.754 
 
  Ускорение на поверхности=   
 
    9.7866945 
 
  Расстояние до центра Земли=   
 
    6382859.7 
 

Неправильно. При чем тут G ? При чем тут радиус Земли?
Вы решили совершенно другую задачу.

В задаче дано  ускорение в двадцати точках.
По этим точкам нужно построить кривую ускорений,  получить уравнение кривой,
по кривой ускорений  построить кривую скоростей, и установить скорость на поверхности Земли.

И интересно: ускорение на поверхности дано:  0         9.8066 (смотри таблицу)
а Вы получили  9.7866945 


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




   

[ssp123]

Неправильно. При чем тут G ? При чем тут радиус Земли?
Вы решили совершенно другую задачу.

В задаче дано  ускорение в двадцати точках.
По этим точкам нужно построить кривую ускорений,  получить уравнение кривой,
по кривой ускорений  построить кривую скоростей, и установить скорость на поверхности Земли.

И интересно: ускорение на поверхности дано:  0         9.8066 (смотри таблицу)
а Вы получили  9.7866945 

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

Давайте подберем функцию, наилучшим образом проходящую через эти точки.
Чтобы не мучаться подгонять полиномом, я предположил ее равной G*M/(h+R)^2

Результаты:
h=6.37637e+06  g=9.80660e+00  расч=9.80713e+00
h=6.37737e+06  g=9.80360e+00  расч=9.80405e+00
h=6.37837e+06  g=9.80050e+00  расч=9.80098e+00
h=6.37937e+06  g=9.79740e+00  расч=9.79790e+00
h=6.38037e+06  g=9.79430e+00  расч=9.79483e+00
h=6.38137e+06  g=9.79120e+00  расч=9.79176e+00
h=6.38237e+06  g=9.78820e+00  расч=9.78870e+00
h=6.38437e+06  g=9.78200e+00  расч=9.78256e+00
h=6.38637e+06  g=9.77590e+00  расч=9.77644e+00
h=6.39137e+06  g=9.76050e+00  расч=9.76115e+00
h=6.39637e+06  g=9.74520e+00  расч=9.74589e+00
h=6.42637e+06  g=9.65420e+00  расч=9.65511e+00
h=6.45637e+06  g=9.56440e+00  расч=9.56559e+00
h=6.47637e+06  g=9.50500e+00  расч=9.50661e+00
h=6.49637e+06  g=9.44700e+00  расч=9.44816e+00
h=6.87637e+06  g=8.45000e+00  расч=8.43277e+00
h=7.37637e+06  g=7.36000e+00  расч=7.32830e+00
h=1.63764e+07  g=1.50000e+00  расч=1.48680e+00
h=5.63764e+07  g=1.25000e-01  расч=1.25457e-01
h=4.06376e+08  g=2.50000e-03  расч=2.41453e-03
Отличаются не сильно? Тогда эту функцию мы имеем полное право загнать в вышеприведенную программу.
Что я и сделал.
То что ускорения в программе не совсем совпадают- возможно , ошибки вычисления. Можно уменьшить dt- будет точнее, и взять h - побольше. Будет точнее, но считать будет компьютер существенно дольше, у меня нетбук, медленный.
То что скорость близка ко второй космической- меня не удивило. Если запустить спутник вертикально вверх, со скоростью чуть меньше второй космической, то он вскоре вернутся, с той же скоростью.

[ssp123]

расчитал еле-еле

Итераций выполнено=   
 
    8793213. 
 
  Скорость на поверхности=                               
 
    11146.854 
 
  Ускорение на поверхности=                             
 
    9.800695 
 
  Расстояние до центра Земли=   
                         
    6378299. 
 

[Вашкевич Виктор]

Давайте подберем функцию, наилучшим образом проходящую через эти точки.
Чтобы не мучаться подгонять полиномом, я предположил ее равной G*M/(h+R)^2

Результаты:
h=6.37637e+06  g=9.80660e+00  расч=9.80713e+00
h=6.37737e+06  g=9.80360e+00  расч=9.80405e+00
h=6.37837e+06  g=9.80050e+00  расч=9.80098e+00
h=6.37937e+06  g=9.79740e+00  расч=9.79790e+00
h=6.38037e+06  g=9.79430e+00  расч=9.79483e+00
h=6.38137e+06  g=9.79120e+00  расч=9.79176e+00
h=6.38237e+06  g=9.78820e+00  расч=9.78870e+00
h=6.38437e+06  g=9.78200e+00  расч=9.78256e+00
h=6.38637e+06  g=9.77590e+00  расч=9.77644e+00
h=6.39137e+06  g=9.76050e+00  расч=9.76115e+00
h=6.39637e+06  g=9.74520e+00  расч=9.74589e+00
h=6.42637e+06  g=9.65420e+00  расч=9.65511e+00
h=6.45637e+06  g=9.56440e+00  расч=9.56559e+00
h=6.47637e+06  g=9.50500e+00  расч=9.50661e+00
h=6.49637e+06  g=9.44700e+00  расч=9.44816e+00
h=6.87637e+06  g=8.45000e+00  расч=8.43277e+00
h=7.37637e+06  g=7.36000e+00  расч=7.32830e+00
h=1.63764e+07  g=1.50000e+00  расч=1.48680e+00
h=5.63764e+07  g=1.25000e-01  расч=1.25457e-01
h=4.06376e+08  g=2.50000e-03  расч=2.41453e-03
Отличаются не сильно? Тогда эту функцию мы имеем полное право загнать в вышеприведенную программу.
Что я и сделал.
То что ускорения в программе не совсем совпадают- возможно , ошибки вычисления. Можно уменьшить dt- будет точнее, и взять h - побольше. Будет точнее, но считать будет компьютер существенно дольше, у меня нетбук, медленный.
То что скорость близка ко второй космической- меня не удивило. Если запустить спутник вертикально вверх, со скоростью чуть меньше второй космической, то он вскоре вернутся, с той же скоростью.

Ясно... Вы решаете совсем другую задачу...
Я лучше с карандашом и на калькуляторе.

[ssp123]

Ясно... Вы решаете совсем другую задачу...
Я лучше с карандашом и на калькуляторе.

Так поставьте её как-то по другому, что ли. Или освойте скайлаб, ваше решение обсудим.

[БЫК]

Ясно... Вы решаете совсем другую задачу...
Я лучше с карандашом и на калькуляторе.

Не надо .
Достаточно чуть-чуть подправить прогу

Код: [Выделить]

x=[0,1000,2000]// вектор высот h (можно донабить точки)
y=[9.8066,9.8036,9.80005] //вектор ускорений
dd=splin(x, y);//вычисляем коэффициенты сплайна
v=0 // Начальная скорость, м/с
dt=1 // Отрезок времени, через который производим вычисления, с
r=400000000 // Расстояние до центра земли, м
cnt=0 // счетчик итераций

while (r>R)
  g=interp(r,x,y,dd);// интерполируем g
  r0=r
  r=r-(v*dt + g*dt*dt*0.5)
  v=v+g*dt
  cnt=cnt+1
end

disp(cnt,"Итераций выполнено=")
disp(v, " Скорость на поверхности= ")
disp(g, " Ускорение на поверхности= ")
disp(r0, " Расстояние до центра Земли= ")

[ssp123]

Приятно, что в скайлабе нашлись сплайны, чем больше про него узнаю, тем больше нравится.
Только вот расчет долго: минут 10 ждал, потом не выдержал и остановил его.
К тому же, результат ничего неожиданного бы не дал.
Вашкевич Виктор, когда переформулируете задачу?

[ssp123]

Что, неужели Вашкевич прав?!


А может, я ошибся?

Код: [Выделить]

stacksize('max');
y=[9.8066, 9.8036, 9.8005, 9.7974, 9.7943, 9.7912, 9.7882, 9.7820, 9.7759, 9.7605, 9.7452, 9.6542, 9.5644, 9.505, 9.447, 8.45, 7.36, 1.50, 0.125, 0.0025];
x=[0, 1000, 2000, 3000, 4000, 5000, 6000, 8000, 10000, 15000, 20000, 50000, 80000, 100000, 120000, 500000, 1000000, 10000000, 50000000, 400000000];
dd=splin(x, y, "monotone");
dr=10000; // Отрезок через который производим вычисления, m
r=0:dr:400000000;
g=interp(r,x,y,dd);// интерполируем g
n=size(r,'c')
for i=1:n
  vv(i)=(2*g(i)*(r(n)-r(i)))^0.5
end
subplot(2,1,1); plot(r,g); xtitle('Ускорение','расстояние от поверхности, м','Ускорение м/с^2')
subplot(2,1,2); plot(r,vv);xtitle('Скорость','расстояние от поверхности, м','Скорость, м/с')

PS. Ошибку нашел, найдет ли Вашеквич свою?

[ssp123]

PPS. Ошибку не нашел... запутался.
Вашкевич, помогайте!

[ssp123]

Ошибку нашел.
Интерполяция сплайнами:

Код: [Выделить]

clear
stacksize('max');
y=[9.8066, 9.8036, 9.8005, 9.7974, 9.7943, 9.7912, 9.7882, 9.7820, 9.7759, 9.7605, 9.7452, 9.6542, 9.5644, 9.505, 9.447, 8.45, 7.36, 1.50, 0.125, 0.0025];
x=[0, 1000, 2000, 3000, 4000, 5000, 6000, 8000, 10000, 15000, 20000, 50000, 80000, 100000, 120000, 500000, 1000000, 10000000, 50000000, 400000000];
dd=splin(x, y, "monotone");
dr=10000; // Отрезок через который производим вычисления, m
r=0:dr:400000000;
g=interp(r,x,y,dd);// интерполируем g
n=size(r,'c')
v(n)=0 // начальная скорость
for i=n-1:-1:1
  v(i)=(2*g(i)*dr+v(i+1)^2)^0.5
end
subplot(2,1,1); plot(r,g); xtitle('Ускорение','расстояние от поверхности, м','Ускорение м/с^2')
subplot(2,1,2); plot(r,v); xtitle('Скорость','расстояние от поверхности, м','Скорость, м/с')

Подгонка к виду const/r^2  и интерполяция

Код: [Выделить]

clear
stacksize('max');
function [zr]=G(c,z)
  zr=(z(2)-c/z(1)^2)^2
endfunction

GG = 6.672d-11
M = 5.976d+24
R = 6.376370d+6

x=[0, 1000, 2000, 3000, 4000, 5000, 6000, 8000, 10000, 15000, 20000, 50000, 80000, 100000, 120000, 500000, 1000000, 10000000, 50000000, 400000000];
y=[9.8066, 9.8036, 9.8005, 9.7974, 9.7943, 9.7912, 9.7882, 9.7820, 9.7759, 9.7605, 9.7452, 9.6542, 9.5644, 9.505, 9.447, 8.45, 7.36, 1.50, 0.125, 0.0025];
x=R+x;
z=[x;y];
c=1.0d+10;

[GM,err]=datafit(G,z,c)
disp(err,"Error=");
dr=10000;
r=R:dr:400000000;
n=size(r,'c');
//g=GM/r^2
for i=1:1:n;
  g(i)=GM/r(i)^2;
end;
v(n)=0 // начальная скорость
for i=n-1:-1:1
  v(i)=(2*g(i)*dr+v(i+1)^2)^0.5
end
subplot(2,1,1); plot(r,g); xtitle('Ускорение','расстояние от поверхности, м','Ускорение м/с^2')
subplot(2,1,2); plot(r,v); xtitle('Скорость','расстояние от поверхности, м','Скорость, м/с')

Похоже, прав Ньютон, а не Вашкевич...

[Вашкевич Виктор]

PPS. Ошибку не нашел... запутался.
Вашкевич, помогайте!

Привет!
Как я понимаю, значения g из таблицы википедии получены по той формуле,
которую задействовали Вы. Так что получается банальность: таблица получена по формуле,
а Вы выбрали формулу для  описания кривой ускорений. Но это ладно, наверное все-таки в
связи с космическими полетами ускорения уточнялись, и тогда можно согласиться с Вашим методом.

Итак, Вы получили скорость на поверхности Земли. А Вы осознаете, что это Вы получили скорость эфира
у поверхности Земли?  Получили ту скорость, которую физики больше ста лет хотят измерить
интерферометром  Майкельсона?


Я то много лет думал, что  измерение скорости эфира интерферометром Майкельсона у физиков
не получается потому, что они исходя из гипотезы "неподвижного эфира"  замеряли скорость
только в горизонтальной плоскости, то есть измеряли скорость поперек потока эфира.
Но  судя по ролику в  теме "эфир есть", последние  замеры делались и по вертикали,
но показаний скорости не получили. Что это значит?
А это значит, что интерферометр не работает.  В общем, получается мировой конфуз:
физики сто лет пользовались неработающим  прибором. Правду сказать, Николай Морозов -
легенда русской науки, - еще в  начале двадцатых годов высказывал предположение,
что интерферометром замерять нельзя, ибо  принцип, заложенный в конструкцию интерферометра  работает.
Ну ладно, с ихним интерферометром.  Но не хотели бы Вы сделать шаги  к знакомству с эфиром?


Скорость втекания эфира в поверхность Земли Вы установили, - это вторая космическая.
И тогда,  если исходить из радиуса Земли  можно узнать площадь поверхности земного шара,
а  зная скорость втекания эфира, легко подсчитать объем втекающего в Землю эфира в одну секунду.
А Ацюковский предлагает  диэлектрическую проницаемость вакуума считать плотностью эфира,
равную 8,85х10^-12 кг. А значит, если объем умножим на плотность, получим  число кг,
- на сколько Земля за счет втекающего эфира увеличивается в одну секунду.

Как? Можете это подсчитать?

[Вашкевич Виктор]

Ошибку нашел.
Интерполяция сплайнами:

Код: [Выделить]

clear
stacksize('max');
y=[9.8066, 9.8036, 9.8005, 9.7974, 9.7943, 9.7912, 9.7882, 9.7820, 9.7759, 9.7605, 9.7452, 9.6542, 9.5644, 9.505, 9.447, 8.45, 7.36, 1.50, 0.125, 0.0025];
x=[0, 1000, 2000, 3000, 4000, 5000, 6000, 8000, 10000, 15000, 20000, 50000, 80000, 100000, 120000, 500000, 1000000, 10000000, 50000000, 400000000];
dd=splin(x, y, "monotone");
dr=10000; // Отрезок через который производим вычисления, m
r=0:dr:400000000;
g=interp(r,x,y,dd);// интерполируем g
n=size(r,'c')
v(n)=0 // начальная скорость
for i=n-1:-1:1
  v(i)=(2*g(i)*dr+v(i+1)^2)^0.5
end
subplot(2,1,1); plot(r,g); xtitle('Ускорение','расстояние от поверхности, м','Ускорение м/с^2')
subplot(2,1,2); plot(r,v); xtitle('Скорость','расстояние от поверхности, м','Скорость, м/с')

Подгонка к виду const/r^2  и интерполяция

Код: [Выделить]

clear
stacksize('max');
function [zr]=G(c,z)
  zr=(z(2)-c/z(1)^2)^2
endfunction

GG = 6.672d-11
M = 5.976d+24
R = 6.376370d+6

x=[0, 1000, 2000, 3000, 4000, 5000, 6000, 8000, 10000, 15000, 20000, 50000, 80000, 100000, 120000, 500000, 1000000, 10000000, 50000000, 400000000];
y=[9.8066, 9.8036, 9.8005, 9.7974, 9.7943, 9.7912, 9.7882, 9.7820, 9.7759, 9.7605, 9.7452, 9.6542, 9.5644, 9.505, 9.447, 8.45, 7.36, 1.50, 0.125, 0.0025];
x=R+x;
z=[x;y];
c=1.0d+10;

[GM,err]=datafit(G,z,c)
disp(err,"Error=");
dr=10000;
r=R:dr:400000000;
n=size(r,'c');
//g=GM/r^2
for i=1:1:n;
  g(i)=GM/r(i)^2;
end;
v(n)=0 // начальная скорость
for i=n-1:-1:1
  v(i)=(2*g(i)*dr+v(i+1)^2)^0.5
end
subplot(2,1,1); plot(r,g); xtitle('Ускорение','расстояние от поверхности, м','Ускорение м/с^2')
subplot(2,1,2); plot(r,v); xtitle('Скорость','расстояние от поверхности, м','Скорость, м/с')

Похоже, прав Ньютон, а не Вашкевич...

А в чем тут, по Вашему, Вашкевич оспаривает Ньютона?

[ssp123]

Здравствуйте, Виктор!

Как я понимаю, значения g из таблицы википедии получены по той формуле,
которую задействовали Вы.

Нет, похоже это экспериментальные значения. Подгоняя их к виду C/r^2 , то есть вычисляя С-константу, я получил ошибку наименьших квадратов, не равную 0.

Итак, Вы получили скорость на поверхности Земли. А Вы осознаете, что это Вы получили скорость эфира
у поверхности Земли?  Получили ту скорость, которую физики больше ста лет хотят измерить
интерферометром  Майкельсона?

Я вообще не могу понять, что такое эфир, и для чего он нужен. Раньше считалось, что он нужен для распространения Э/М волн. Я легко представляю их распространение в вакууме, точно так же как распространение материи в нем.
То есть, не осознаю. Но допускаю его существование, пока не доказано обратное.

Но не хотели бы Вы сделать шаги  к знакомству с эфиром?

С удовольствием, благо пока нахожусь в ЧАО, и времени свободного много.

Как? Можете это подсчитать?

Легко.

А в чем тут, по Вашему, Вашкевич оспаривает Ньютона?

Вы вычисляли вторую космическую формулу, и у Вас она , помнится, получалась на порядок больше классической.
Значит, у Вас была ошибка.
А Семенов (это я- Семенов, ssp) - с помощью компьютера, тремя различными способами, практически подтвердил старика Ньютона.
 

[ssp123]

Масса втекаемого за 1 секунду в Землю эфира
m=4*%pi*(6.376370d+6)^2*8.85d-12*11200
 m  = 50642855.  кг