Дискретно-аналитический метод решения задачи.
Для решения задачи будем использовать дискретно-аналитический метод, который состоит в следующем: по оси x осуществляется конечно-разностная аппроксимация, а по оси времени t рассматривается непрерывная задача (рис. 4.4).
Рис. 4.4. Схема дискретизации.
Введем обозначения:
; 
, 
, (4.55)
где в простейшем случае
. (4.56)
Здесь 
– количество внутренних узлов конечно-разностной сетки, причем пусть – нечетное число.
Для всех внутренних узлов 
получим конечно-разностное уравнение – дискретный аналог уравнения колебаний (4.54):
, 
, (4.57)
где
(4.58)
– вторая конечная разность, приближенно представляющая вторую производную от искомой функции 
по аргументу 
.
В соответствии с краевыми условиями из (4.54) для граничных узлов, очевидно, можем записать:
; 
; 
; 
. (4.59)
C учетом (4.59) преобразуются уравнения (4.57) Имеем:
; (4.60)
; (4.61)
, 
; (4.62)
; (4.63)
. (4.64)
Обоснуем, например, (4.60) и (4.61):
;
Введя обозначение
, (4.65)
можем представить разрешающую систему конечно-разностных уравнений (4.60)-(4.64) в матричном виде
(4.66)
где
. (4.67)
Заметим, что матрица 
положительно определена, т.е. все ее собственные числа положительные (в этом можно убедиться при их непосредственном вычислении).
Общее решение задачи (4.66) имеет вид:
. (4.68)
По условию рассматриваемой задачи
, (4.69)
откуда следует, что
, где 
. (4.70)
Это значит, что в векторе 
лишь один «срединный элемент» (с номером 
) равен единице, а остальные элементы равны нулю. Ненулевой элемент вектора соответствует узлу конечно-разностной сетки с координатой 
, в котором в момент времени 
приложено сосредоточенное ударное воздействие величиной 
.
Подставив (4.69) в (4.68), получим окончательный вид общего решения:
. (4.71)
Варианты задания.
– величина приложенного сосредоточенного ударного воздействия; 
; 
; 
– номер группы, 
– номер студента по журналу.
Принять количество количество внутренних узлов конечно-разностной сетки 
.
Пример соответствующего M-файла (ниже задано 
, 
):
| format bank n=input('введите n='); G=3; S=12; L=300; P=300; h=L/(n+1); alfa=10^8*(100+G+S); x=0:h:L; a0=6*eye(n); a0(1,1)=5; a0(n,n)=5; a1=ones(n-1,1); a2=ones(n-2,1); A=a0-4*(diag(a1,-1)+diag(a1,1))+diag(a2,-2)+diag(a2,2) A=alfa*A/h^4; [T,J]=eig(A) T_INV=inv(T) for i=1:n fJ(i,i)=sqrt(J(i,i)); end sq_A=T*fJ*T'; t0=pi/(4*fJ(n,n)); tmax=125*t0; kmax=7; ic=(n+1)/2; tk=[t0,2*t0,3*t0,tmax/2,tmax-2*t0,tmax-t0,tmax]; res=zeros(kmax,n+2); for k=1:kmax t=tk(k); for i=1:n fJt(i,i)=sin(fJ(i,i)*t); end Y_t=-P*inv(sq_A)*T*fJt*T'; res(k,2:n+1)=Y_t(:,ic)'; end result_y=[tk' res] hold on plot(x,res(1,1:n+2),'r') plot(x,res((kmax+1)/2,1:n+2)) plot(x,res(kmax,1:n+2),'g') grid on |
Результаты расчета:
A =
5.00 -4.00 1.00 0 0 0 0
-4.00 6.00 -4.00 1.00 0 0 0
1.00 -4.00 6.00 -4.00 1.00 0 0
0 1.00 -4.00 6.00 -4.00 1.00 0
0 0 1.00 -4.00 6.00 -4.00 1.00
0 0 0 1.00 -4.00 6.00 -4.00
0 0 0 0 1.00 -4.00 5.00
T =
0.19 0.35 0.46 0.50 0.46 -0.35 0.19
0.35 0.50 0.35 0.00 -0.35 0.50 -0.35
0.46 0.35 -0.19 -0.50 -0.19 -0.35 0.46
0.50 -0.00 -0.50 -0.00 0.50 0.00 -0.50
0.46 -0.35 -0.19 0.50 -0.19 0.35 0.46
0.35 -0.50 0.35 0.00 -0.35 -0.50 -0.35
0.19 -0.35 0.46 -0.50 0.46 0.35 0.19
J =
187.52 0 0 0 0 0 0
0 776.35 0 0 0 0 0
0 0 12333.06 0 0 0 0
0 0 0 32363.46 0 0 0
0 0 0 0 61872.89 0 0
0 0 0 0 0 94314.02 0
0 0 0 0 0 0 119787.27
T_INV =
0.19 0.35 0.46 0.50 0.46 0.35 0.19
0.35 0.50 0.35 -0.00 -0.35 -0.50 -0.35
0.46 0.35 -0.19 -0.50 -0.19 0.35 0.46
0.50 0.00 -0.50 -0.00 0.50 0.00 -0.50
0.46 -0.35 -0.19 0.50 -0.19 -0.35 0.46
-0.35 0.50 -0.35 0.00 0.35 -0.50 0.35
0.19 -0.35 0.46 -0.50 0.46 -0.35 0.19
result_y =
0.00 0 0.00 0.00 -0.02 -0.65 -0.02 0.00 0.00 0
0.00 0 0.00 0.03 -0.13 -1.16 -0.13 0.03 0.00 0
0.01 0 0.02 0.07 -0.39 -1.43 -0.39 0.07 0.02 0
0.14 0 -1.87 -3.63 -4.97 -4.68 -4.97 -3.63 -1.87 0
0.28 0 1.04 2.20 3.47 3.47 3.47 2.20 1.04 0
0.28 0 1.05 2.61 3.44 3.46 3.44 2.61 1.05 0
0.28 0 1.14 2.98 3.37 3.50 3.37 2.98 1.14 0
Рис. 4.5. Колебания балки при 
, 
и 
.
Здесь, как и в предыдущей лабораторной работе, 
– порядковый номер точки табуляции решения по времени 
.
Дата добавления: 2015-10-09; просмотров: 853;