Решение дифуров методом Рунге-Кутта

[Toolen]

Yu_K, мне не нужны солитоны, мне нужно токько получить три графика зависимости координат грузов, их кинетической и потенциальной энергии от времени

[Yu_K]

Цитата:

Есть три дифура в рамка физической задачи, их необходимо решить

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

А солитоны - это так для кругозора - особо не предлаем - вот здесь есть ссылка на оригинальную работу ФПУ

http://www.sciyouth.ru/sites/11/Researc ... eferat.htm (я было думал в Самаре всплеск интереса к солитонам ).

[ewert]

Yu_K в сообщении #162520 писал(а):

Можно решать аналитически - трехдиагональная матрица - и как заметил ewert искать ее собственные числа и затем выписать решение в виде линейной комбинации гармоник.

Ну, буквально-то этого я не замечал. Там матрица 6 на 6; её диагональные блоки (при соответствующей записи) -- нулевые, один из внедиагональных -- и впрямь трёхдиагонален и симметричен, другой представляет собой единичную матрицу.

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

[Yu_K]

Там ошибка у Toolen в записи уравнений - матрица трехдиагональная - см. например Пановко ЯГ - Введение в теорию колебаний,1991, стр 80 уравнение (4.25).

[ewert]

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

[Toolen]

Yu_K, это моя работа

[Yu_K]

Я так подумал - случайное совпадение четырех параметров - музыка, Самара, Ферми Паста Улам, университет - слишком маловероятно. А с Рунге-Куттом проблемы кончились? Каким-нибудь пакетом владеете математическим? Тут в принципе есть одна проблемка - как автоматизировать постановку задачи для ПК для произвольного кол-ва шариков N, чтобы не набирать каждый раз при каждом N новое кол-во уравнений. В Маткаде это решается достаточно просто.

[Toolen]

Если бы я был математиком, то не полез бы на этот форум, сейчас я музыкант (с дипломом АСУшника), но учеба заставляет время от времени программировать всякую всячину. Есть код на МатЛабе, но почему-то не рабочий, а в синтаксисе МатЛаба я не разбераюсь, могу перслать код.

[Yu_K]

Начните, например, с двух шариков и три пружины- будет 3 уравнения второго порядка - от них перейдите к системе первого порядка - будет 6 уравнений - а дальше можно любой мат пакет - там везде есть примеры решения задачи Коши для систем ОДУ - и как правило есть встроенный Рунге-Кутт. Хотя нет проблем метод запрограммировать самому.

[Toolen]

Это и нужно сделать, но для начала хотелось бы разобраться напримере.

Сопсно код для МАТЛАБа:

function fpu1
function dy=fpu1(t,y);
N=32;
alpha=0.25;
D(N+1)=y(2)-2*y(1)+alpha*((y(2)-y(1))\verb^2-y(1)\verb^2);
D(1)=y(N+1);
D(2*N)=y(N-1)-2*y(N)+alpha*(y(N)\verb^2-(y(N)-y(N-1))\verb^2);
D(N)=y(2*N);
for I=2:N-1,
D(N+I)=y(I+1)+y(I-1)-2*y(I)+alpha*((y(I+1)-y(I))\verb^2-(y(I)-y(I-1))\verb^2);
D(I)=y(N+I);
end
dy=D';
N=32; %Number of particles must be a power of 2
alpha=0.25; %Nonlinear parameter
TMAX=10000;
DT=20;
tspan=[0:DT:TMAX];
options=odeset('Reltol',1e-4,'OutputFcn','odeplot','OutputSel',[1,2,N]); % Test different tolerances, changing Reltol
for I=1:N,
a=1; b(I)=a*sin(pi*I/(N+1)); b(I+N)=0; % FPU initial condition
%a=1; b(I)=a*sin(pi*N*I/(N+1)); b(I+N)=0; % Zabusky-Deem init. cond.
%k=0.8; sk=(sinh(k))\verb!^!2; ek=exp(-k); i1=I-N/4; i2=i1-N/2; %Solitons
%b(I)=-0.5/alpha*log((1+exp(2*k*(i1-1)))/(1+exp(2*k*i1))); % Kink
%b(I)=b(I)+0.5/alpha*log((1+exp(2*k*(i2-1)))/(1+exp(2*k*i2))); % Anti-kink
%b(I+N)= sk*ek/alpha/cosh(k*i1)/(exp(-k*i1)+exp(k*i1)*exp(-2*k));
%b(I+N)=b(I+N)-sk*ek/alpha/cosh(k*i2)/(exp(-k*i2)+exp(k*i2)*exp(-2*k));
omegak2(I)=4*(sin(pi*I/2/N))\verb^2; % Mode Frequencies
end
[T,Y]=ode45('fpu1',tspan,b',options,N); % Time integration
for IT=1:(TMAX/DT),
TIME(IT)=IT*DT*sqrt(omegak2(1))/2/pi; % Time iteration loop
YX(IT,1:N+1)=[0 Y(IT,1:N)];
YV(IT,1:N+1)=[0 Y(IT,N+1:2*N )];
sXF(IT,:)=imag(fft([YX(IT,1:N+1) 0 -YX(IT,N+1:-1:2)]))/sqrt(2*(N+1));
sVF(IT,:)=imag(fft([YV(IT,1:N+1) 0 -YV(IT,N+1:-1:2)]))/sqrt(2*(N+1));
Energ(IT,1:N)=(omegak2(1:N).*(sXF(IT,2:N+1).\verb^2)+sVF(IT,2:N+1).\verb^2)/2;
for J=2:N-1, % Space loop
DifY(IT,J)=Y(IT,J+1)-Y(IT,J);
end
end
plot(TIME,Energ(:,1),TIME,Energ(:,2),TIME,Energ(:,3),TIME,Energ(:,4));
surf(DifY); % Space derivative field to show the soliton dynamics

[TOTAL]

Toolen писал(а):

Это и нужно сделать, но для начала хотелось бы разобраться напримере.

С чем разобраться?
Метод Вам рассказали. Программируйте сами.