4.6.4. Алгоритм решения краевой задачи методом прогонки.

1. Выбирается h, и находятся узловые точки X_i;

2. В узловых точках находятся коэффициенты М_i, К_i и F_i (при необходимости);

3. По (4.68) находятся прогоночные коэффициенты в начальной точке: С₀ и D₀;

4. По (4.70) рассчитываются прогоночные коэффициенты C_i и D_i;

5. По (4.72) находится значение Y_n;

6. Используя формулы (4.73) или (4.66), последовательно вычисляют значения

Y_n-1, Y_n-2 ,.., Y₁.

7. Значение Y₀ находится по (4.74).

8. Для проверки точности решение величина шага изменяется вдвое, и алгоритм повторяется, начиная с п.1;

9. Оценивается точность найденного решения сравнением значений Y_i и Y_i. Если найденная погрешность удовлетворяет заданной точности, решение завершается. В противном случае вычисления повторяются с уменьшенным вдвое шагом.

Лабораторная работа 5. Определение интеграла функции.

Задание:

1. Вычислить определенный интеграл (n – количество отрезков деления участка интегрирования):

n=10

n=4

n=10

n=11

Вычисления каждого интеграла производить тремя способами:

1. По формуле Ньютона-Лейбница:

где y₁,…,у_n – значения функции f(x) в точках х₁,…,х_n;

h=(b-a)/n; [a,b] – участок интегрирования, n – количество отрезков деления участка интегрирования.

2. По формуле трапеций:

3. По формуле Симпсона:

, где n=2m

2. Вычислить интеграл с точностью е=0,0001:

методом трапеций

методом прямоугольников

методом трапеций

методом трапеций

Вычисления выполнять в одном цикле, при достижении нужной точности для одного интеграла вычисления для другого продолжать.

3. Составить программу вычисления определенного интеграла с заданной точностью:

Вари-ант зада-ния	Подынте-гральная функция	Метод численного решения	Число отрезков n	Интервал интегри-рования [a,b]	Требуе-мая точность 
1	ln^2x/x	Трапеций	60	[1;4]	10^-4
2		Прямоугольни-ков	50	[1;2,5]	0,5*10^-3
3	x^x(1+lnx)	Трапеций	40	[1;3]	10^-4
4	Cos x	Трапеций	60		10^-4
5	Sin^2x	Трапеций	60		0,5*10^-3
6	xe^xsinx	Трапеций	100	[0;1]	10^-4
7	(ln x/x)^2	Прямоугольни-ков	50	[1;2,5]	10^-4
8	x arctg x	Трапеций	50	[0;3]	0,5*10^-3
9		Прямоугольни-ков	100	[0;2]	10^-5
10	e^xcos^2x	Трапеций	60	[0; ]	10^-4
11	x^3/(3+x)	Прямоугольни-ков	80	[1;2]	0,5*10^-4
12	(ln x/x)^3	Трапеций	50	[1;2]	10^-4
13		Прямоугольни-ков	50	[0;2]	10^-4
14	x^2sin2x	Трапеций	100	[1;2]	10^-4
15		Трапеций	50	[1;2]	0,5*10^-3

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

Построить зависимость разности значений интеграла для двух последовательных чисел отрезков разбиения от числа отрезков разбиения.

Пример программы для варианта 1. Исходные данные: число интервалов разбиения участка интегрирования –ячейка «Е4», интервал интегрирования – ячейки «F4» и «F5», заданная точность вычисления – ячейка «G4». Результаты расчета – в столбце «D».

Sub Integral()

N = Range("E4")

Xmin = Range("F4")

Xmax = Range("F5")

eps = Range("G4")

i = 1

Integro = 0

dX = (Xmax - Xmin) / N

Do

Integro2 = Integro

Integr = (Xmin + Xmax) / 2

For X = Xmin + dX To Xmax - dX Step dX

Integro = Integro + F(X)

Next X

Integro = dX * Integro

Range("D" + LTrim$(Str$(10 + i))) = Integro

Range("E" + LTrim$(Str$(10 + i))) = Integro2

Range("F" + LTrim$(Str$(10 + i))) = Abs(Integro2 - Integro)

i = i + 1

dX = dX / 2

Loop While Abs(Integro - Integro2) > eps

End Sub

Function F(X)

F = Log(X) ^ 2 / X

End Function

Подпрограмма-функция используется для вычисления значения исследуемой функции, имеет один параметр – значение независимой переменной. Для вызова подпрограммы-функции достаточно в вызывающей программе указать имя подпрограммы-функции с ее параметрами.

Лабораторная работа 6. Решение нелинейных уравнений.

Задание:

1. Решить методом хорд уравнение:

x³-3x²+x-2=0

с точностью до 10^-4, используя формулу

x₁=a-(b-a)f(a)/(f(b)-f(a)),

где f(a) и f(b) – значения функции на концах интервала [a,b], в котором лежит уточняемый корень уравнения.

Кроме того, функция f(х) непрерывна на отрезке [a,b] вместе со своими производными первого и второго порядков и производные сохраняют знак и не обращаются в 0 на всем участке. При решении уравнения методом хорд или касательных нужно определить знаки функции в различных точках и найти интервал [a,b] изменения знака функции. В качестве начального интервала выделения корня следует взять интервал [2,3].

2. Решить методом касательных уравнение:

x³-2x²+x-3=0

с точностью до 10^-4, используя формулу

x=b-f(b)/f`’(b),

или x=a-f(a)/f`’(a).

Выбор точки a или b производится, исходя из условия равенства знака функции и знака второй производной в данной точке.

3. Решить методом хорд и касательных уравнение

(4+x)²(е^х-е^-х)=18.

В качестве начального интервала выделения корня следует взять интервал [1,2;1,3].

Пример программы для варианта 1. Исходные данные: начальный интервал выделения корня – ячейки «С8» и «С9», заданная точность вычисления – ячейка «С10». Результаты расчета корня – в столбце «E», модуль значения функции - в столбце «F».

Sub Метод_хорд()

a = Range("C8")

b = Range("C9")

eps = Range("C10")

X1 = a

i = 1

Do

X0 = X1

X1 = a - (b - a) * F(a) / (F(b) - F(a))

If (F(X1) > 0 And F(a) > 0) Or (F(X1) < 0 And F(a) < 0) Then

a = X1

Else

b = X1

End If

Range("E" + LTrim$(Str$(7 + i))) = X1

Range("F" + LTrim$(Str$(7 + i))) = F(X1)

i = i + 1

Loop While Abs(X1 – X0) > eps

End Sub

Function F(X)

F = X ^ 3 - 3 * X ^ 2 + X - 2

End Function

Пример программы для варианта 1. Решение методом касательных. Используются три подпрограммы-функции: 1 – вычисление значения функции по значению переменной, 2 - вычисление значения производной функции в точке при принятом приращении переменной dX, 3 – определение знака второй производной функции.

Sub Метод_касательных()

a = Range("C8")

b = Range("C9")

eps = Range("C10")

dX = Range("C11")

i = 1

Do

If (F(b) > 0 And ddF(b, dX) > 0) Or (F(b) < 0 And ddF(b, dX) < 0) Then

DX = F(b) / dF(b, dX)

X = b - DX

b = X

Else

DX = F(a) / dF(a, dX)

X = a - DX

a = X

End If

Range("E" + LTrim$(Str$(7 + i))) = X

Range("F" + LTrim$(Str$(7 + i))) = F(X)

i = i + 1

While Abs(DX) > eps

End Sub

Function F(X)

F = X ^ 3 - 3 * X ^ 2 + X - 2

End Function

Function dF(X, dX)

F1 = X ^ 3 - 3 * X ^ 2 + X - 2

F2 = (X + dX) ^ 3 - 3 * (X + dX) ^ 2 + (X + dX) - 2

dF = (F2 - F1) / dX

End Function

Function ddF(X, dX)

F1 = X ^ 3 - 3 * X ^ 2 + X - 2

F2 = (X + dX) ^ 3 - 3 * (X + dX) ^ 2 + (X + dX) - 2

F3 = (X - dX) ^ 3 - 3 * (X - dX) ^ 2 + (X - dX) - 2

ddF = Abs(F2 - F1) - Abs(F1 - F3)

End Function

Лабораторная работа 7. Решение систем линейных уравнений.

Задание:

Решить методом Гаусса систему уравнений:

1		2х1+3х2+11х3+5х4=2; х1+х2+5х3+2х4=1; 2х1+х2+3х3+2х4=-3; х1+х2+3х3+4х4=-3.
2		2х1+х2-0,1х3+х4=2,7; 0,4х1+0,5х2+4х3-8,5х4=21,9; 0,3х1-х2+х3+5,2х4=-3,9; х1+0,2х2+2,5х3-х4=9,9.
3		0,15х1+2,11х2+30,75х3=-26,38; 0,64х1+1,21х2+2,05х3=1,01; 3,21х1+1,53х2+1,04х3=5,23.
4		1,15х1+0,42х2+100,71х3=-198,70; 1,19х1+0,55х2+0,82х3=2,29; 1,00х1+0,35х2+3,00х3=-0,65.
5		2,6х1-4,5х2-2,0х3=19,07; 3,0х1+3,0х2+4,3х3=3,21; 6,0х1+3,5х2+3,0х3=-18,25.

Методические указания. В основе метода Гаусса лежит идея последовательного исключения неизвестных. Для примера рассмотрим систему четырех уравнений с четырьмя неизвестными: .

а11х1+а12х2+а13х3+а14х4=а15; а21х1+а22х2+а23х3+а24х4=а25; а31х1+а32х2+а33х3+а34х4=а35; а41х1+а42х2+а43х3+а44х4=а45.

Пусть а₁₁0. Разделив первое уравнение системы на а₁₁, получим

х₁+b₁₂х₂+b₁₃х₃+b₁₄х₄=b₁₅;

где b_1j=a_1j/a₁₁; (j=2,3,4,5)

Используя это уравнение, исключим х₁ из второго, третьего, четвертого уравнений системы. Для этого необходимо умножить это уравнение на а₂₁, а₃₁, а₄₁ и вычесть соответственно из второго, третьего и четвертого уравнения системы.

После выполнения действий получим:

а22^1х2+а23^1х3+а24^1х4=а25^1; а32^1х2+а33^1х3+а34^1х4=а35^1; а42^1х2+а43^1х3+а44^1х4=а45^1.

где а_ij¹=а_ij-а_i1b_1j (i=2,3,4; j=2,3,4,5).

Затем первое уравнение полученной системы делим на а₂₂¹ и получим

х₂+b₂₃¹х₃+b₂₄¹х₄=b₂₅¹;

где b_2j¹=a_2j¹/a₂₂¹; (j=3,4,5)

Исключив х₂ из системы уравнений, получим:

а33^2х3+а34^2х4=а35^2; а43^2х3+а44^2х4=а45^2.

где а_ij²=а_ij¹-а_i2¹b_2j¹ (i=3,4; j=3,4,5).

Разделив первое уравнение последней системы на а₃₃², получим

х₃+b₃₄²х₄=b₃₅²;

где b_3j²=a_3j²/a₃₂² ; (j=4,5)

Затем исключаем х₃ из второго уравнения системы и получаем

а₄₄³х₄=а₄₅³

где а_4j³=а_4j²-а₄₃²b_3j² (j=4,5).

После всех преобразований получим систему треугольного вида:

x1+b12х2+b13х3+ b14х4 = b15; х2+b23^1х3+b24^1х4=b25^1; х3+b34^2х4=b35^2; а44^3х4=а45^3.

откуда последовательно получаем:

х₄=а₄₅³/а₄₄³;

х₃=b₃₅²-b₃₄²х₄;

х₂=b₂₅¹-b₂₃¹х₃-b₂₄¹х₄;

x₁=b₁₅-b₁₂х₂-b₁₃х₃-b₁₄х_4.

Процесс нахождения решения по методу Гаусса распадается на два этапа. Первый этап — приведение системы к треугольному виду (прямой ход). Второй этап— вычисление неизвестных по полученным формулам принято называть обратным ходом.

Пример программы решения системы линейных уравнений методом Гаусса.

Количество уравнений – до 4. Коэффициенты уравнений располагаются, начиная с ячейки С2. По количеству коэффициентов в столбце «С» программа определяет число уравнений и столбцы расположения следующих коэффициентов. Результаты нахождения неизвестных располагаются в строке, начиная с ячейки «С7» для неизвестной при коэффициентах, расположенных в столбце «С».

Sub Линейная_система2()

Dim A(4, 5) As Double

Dim B(4, 5) As Double

Dim X(4) As Double

Dim Buc$(5)

Buc$(0) = "C": Buc$(1) = "D": Buc$(2) = "E":

Buc$(3) = "F": Buc$(4) = "G"

'

'

'Чтение коэффициентов уравнений

i = 1

While Range("C" + LTrim$(Str$(i + 1)))

' Коэффициенты 1-го столбца

A(i - 1, 0) = Range("C" + LTrim$(Str$(i + 1)))

i = i + 1

Wend

N = i - 1 'Число уравнений (неизвестных)

For j = 1 To N

For k = 1 To N

A(j - 1, k) = Range(Buc$(k) + LTrim$(Str$(j + 1)))

Next k

Next j

'

'

' Последовательное исключение неизвестных

For k = 1 To N - 1

For j = k To N

B(k - 1, j) = A(k - 1, j) / A(k - 1, k - 1)

Next j

For j = k To N

For i = k To N - 1

A(i, j) = A(i, j) - A(i, k - 1) * B(k - 1, j)

Next i

Next j

Next k

'

'

' Нахождение значений неизвестных

X(N - 1) = A(N - 1, N) / A(N - 1, N - 1)

For i = N - 2 To 0 Step -1

S = 0

For j = i + 1 To N - 1

S = S + B(i, j) * X(j)

Next j

X(i) = B(i, N) - S

Next i

'

' Печать значений неизвестных

For k = 1 To N

Range(Buc$(k - 1) + "7") = X(k - 1)

Next k

End Sub

Для чтения из ячеек Листа и записи в ячейки двумерных массивов часто удобнее пользоваться функцией Cells(i,j), позволяющей обращаться к ячейке, расположенной на i-й строке и j-м столбце. Напримр, значение из ячейки В3 может быть прочитано оператором:

Alfa = Cells(3,2)

Макрос для чтения таблицы данных, расположенных, начиная с ячейки А1 может иметь вид:

Sub WriteDate()

'

Dim A(10, 10) As Single

Set Data = ActiveSheet.UsedRange

For i = 1 To Data.Rows.Count

For j = 1 To Data.Columns.Count

A(i - 1, j - 1) = Cells(i, j)

Next j

Next i

'

End Sub

Лабораторная работа 8. Решение обыкновенного дифференциального уравнения методом Рунге-Кутта.

Обыкновенное дифференциальное уравнение – уравнение, содержащее производные неизвестной функции от одной независимой переменной.

Задание:

1) Найти решение уравнения:

у=0,25у²+х²

с начальными условиями

у(0)=-1; h=0,1 на отрезке [0;0,5].

2) у=x+y^0,5

с начальными условиями

у(0,5)=0,724; h=0,1 на отрезке [0,5;1,5].

3) у=0,5xy

с начальными условиями

у(0)=1; h=0,1 на отрезке [0;1].

4) у=2x+cosy

с начальными условиями

у(0)=0; h=0,02 на отрезке [0;0,1].

5) у=y²+2/x

с начальными условиями

у(0)=1; h=0,2 на отрезке [1;1,7].

6) у=y³+x³

с начальными условиями

у(0)=0,5; h=0,1 на отрезке [0;0,8].

Методические указания. Обыкновенным дифференциальным уравнением называется уравнение, содержащее производные неизвестной функции от одной независимой переменной.

Наивысший порядок производной определяет порядок дифференциального уравнения. Дифференциальное уравнение n-й степени можно свести к n уравнениям первой степени. Существует ряд аналитических методов решения дифференциальных уравнений, но часто бывают случаи, когда трудно, а иногда и невозможно их применить. В таких ситуациях используются численные методы решения дифференциальных уравнений с использованием ЭВМ. Существует много методов численного интегрирования дифференциальных уравнений, но мы остановимся на методе Рунге— Кутта, в котором значение функции y_i+1 в точке x_i+1 вычисляется при помощи значения функции у_i в точке x_i.

Рассмотрим задачу Коши для дифференциального уравнения

y=f(x,y)

с начальными условиями y(x₀)=y₀.

Обозначим через y_i приближенное значение решения в точке x_i. По методу Рунге-Кутта вычисляют приближенные значения y_i+1 в точке

x_i+1=x_i+h, где h – шаг.

Вычисляют по формулам:

y_i+1=y_i+y_i;

y_i=(k₁ⁱ+2 k₂ⁱ+2k₃ⁱ+ k₄ⁱ)/6,

где

k₁ⁱ=hf(x_i;y_i);

k₂ⁱ=hf(x_i+h/2;y_i+k₁ⁱ/2);

k₃ⁱ=hf(x_i+h/2;y_i+k₂ⁱ/2);

k₄ⁱ=hf(x_i+h;y_i+k₃ⁱ);

Все вычисления можно расположить в таблице:

I	x	Y	k=hf(x,y)	y
0	x0	y0	k1^0	k1^0
	x0+h/2	y0+k1^0/2	k2^0	2k2^0
	x0+h/2	y0+k2^0/2	k3^0	2k3^0
	x0+h	y0+k3^0	k4^0	k4^0
				y0
1	x1	y1

Порядок заполнения таблицы следующий.

1. В первой строке таблицы записывают данные значения x₀ и y₀.

2. Вычисляют f(x_o,у_о), умножают на h и заносят в таблицу в виде k₁⁰.

3. Во второй строке таблицы вычисляют и записывают х_о+h/2 и y_o+k₁⁰/2.

4. Вычисляют значение f(x₀+h/2, y₀+k₁⁰/2) и после умножения на h записывают в таблицу в качестве k₂⁰.

5. В третьей строке записывают х_о+h/2 и y_o+k₂⁰/2.

6. Вычисляют значение f(x₀+h/2, y₀+k₂⁰/2) и после умножения на h записывают в таблицу в виде k₃⁰.

7. В четвертой строке записывают х_о+h и y_o+k₃⁰.

8. Вычисляют f(x₀+h, y₀+k₃⁰) и после умножения на h заносят в таблицу в виде k₄⁰.

9. В столбец y записывают числа k₁⁰, 2k₂⁰, 2k₃⁰, k₄⁰.

10. Суммируют числа в столбце y, результат делят на 6 и заносят в качестве y₀.

11. Вычисляют у₁=у₀+у₀.

Затем продолжают в том же порядке, но за начальную точку берут (х₁,у₁).

Пример программы решения линейного дифференциального уравнения первого порядка вида dT/dt=(Tg-T)/eps (уравнение динамики одноемкостного звена).

Исходные данные берутся из ячеек активного листа D7:D12, результаты расчета размещаются в столбцы F, G, H.

Sub Инерц_звено()

tau0 = Range("D7")

T0 = Range("D8")

eps = Range("D9")

Tg = Range("D10")

dtau = Range("D11")

tau = Range("D12")

T = T0

i = 1

For TimeProcess = tau0 To tau Step dtau

k1 = F(Tg + 10 * i, T, eps) * dtau

k2 = F(Tg + 10 * i, T + k1 / 2, eps) * dtau

k3 = F(Tg + 10 * i, T + k2 / 2, eps) * dtau

k4 = F(Tg + 10 * i, T + k3, eps) * dtau

dT = (k1 + 2 * k2 + 2 * k3 + k4) / 6

T = T + dT

i = i + 1

Range("F" + LTrim$(Str$(i + 6))) = TimeProcess

Range("G" + LTrim$(Str$(i + 6))) = T

Range("H" + LTrim$(Str$(i + 6))) = Tg + 10 * (i - 1)

Next TimeProcess

End Sub

Function F(Tg, T, eps)

F = (Tg - T) / eps

End Function

В программе используется подпрограмма-функция.

Лабораторная работа 9. Построение эмпирической формулы методом наименьших квадратов.

Задание:

Для набора пар данных (независимая переменная Х, функция У) методом наименьших квадратов построить линию регрессии:

Х	1	2	3	4	5	6	7	8	9	10
У	6	4	12	7	22	14	23	16	10	20

1. используя в качестве линии регрессии прямую линию.

2. используя в качестве линии регрессии квадратичную параболу.

Пример программы линейного регрессионного анализа.

Исходные данные берутся из столбцов активного листа С и D, результаты расчета размещаются в столбец Е, коэффициенты регрессии – в ячейки В1 и В2.

Sub Линейная_регрессия()

Dim X(200), Y(200)

i = 1

While Range("C" + LTrim$(Str$(i)))

X(i) = Range("C" + LTrim$(Str$(i)))

Y(i) = Range("D" + LTrim$(Str$(i)))

i = i + 1

Wend

Limit = i - 1

SummaX = 0

SummaY = 0

SummaX2 = 0

SummaXY = 0

For j = 1 To Limit

SummaX = SummaX + X(j)

SummaY = SummaY + Y(j)

SummaX2 = SummaX2 + X(j) ^ 2

SummaXY = SummaXY + X(j) * Y(j)

Next j

a1 = (Limit * SummaXY - SummaY * SummaX) / (Limit * SummaX2 - SummaX ^ 2)

a0 = (SummaY - a1 * SummaX) / Limit

Range("B1") = a0

Range("B2") = a1

For k = 1 To Limit

Range("E" + LTrim$(Str$(k))) = a0 + a1 * X(k)

Next k

End Sub

Лабораторная работа 10. Моделирование кривошипно-ползунного механизма

Рис. 10.1. Схема кривошипно-ползунного механизма

В приводимой ниже задаче анализа движения кривошипно-ползунного механизма (рис.10.1) требуется определить параметры движения: координаты и перемещения звеньев, скорости, угловые скорости и ускорения звеньев, направления скоростей для приведения сил, используя моделирование последовательных положений механизма.

Алгоритм моделирования движения кривошипно-ползунного механизма включает в себя следующие этапы.

Расчет характеристик положения. Пусть угол ₁ поворота ведущего звена связан с временем t так: ₁=₁t, где ₁ — угловая скорость ведущего звена. Задавая ₁=2n с^-1, установим значение шага t в зависимости от приращений угла ₁: t=₁/₁. Полагая направление осей положительным, как показано на рис. 10.1, положение точки В определим ее координатами:

Положение точки С — центра ползуна — определим, зная длину шатуна l_BC из системы уравнений:

откуда

Координаты точки S, расположенной на звене ВС на расстоянии l_BS от точки В, определяются из подобных треугольников СВВ' и CSS':

Для определения угла  поворота шатуна с целью дальнейшего определения скоростей и ускорений целесообразно вычислить =arctg y_S/(x_C-x_S). Все эти величины следует сохранять в соответствующих массивах, если в дальнейшем планируется применять процедуры дифференцирования в интегрирования (например, стандартные).

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

Для получения скоростей V_B, V_C, V_S в выделенных точках В, С, S достаточно продифференцировать координаты х_i, у_i, т.е. получить v_xi=dx_i/dt; v_yi=dy_i/dt.

Чтобы получить полные скорости и их направления, необходимо для i-го звена механизма во всех j-x точках положения механизма вычислить

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

Расчет приведенных сил. Для расчета приведенных сил каждой j-й точки k-то звена механизма пользуемся формулами

Пример программы расчета координат и линейных скоростей точек В, С, и S кривошипного механизма.

Sub Ползун()

Lab = Range("C14")

Lbc = Range("C15")

Lbs = Range("C16")

Omega1 = Range("C17")

dFi1 = Range("C18")

Lcs = Lbc - Lbs

Pi = 3.14

dFi1 = dFi1 * Pi / 180

dTau = 0.0001

Period = 2 * Pi / Omega1

i = 2

For Fi1 = 0 To 2 * Pi Step dFi1

Tau = Fi1 / Omega1

Xb = FXb(Lab, Omega1, Tau)

Xb1 = FXb(Lab, Omega1, Tau + dTau)

VXb = (Xb1 - Xb) / dTau

Yb = FYb(Lab, Omega1, Tau)

Yb1 = FYb(Lab, Omega1, Tau + dTau)

VYb = (Yb1 - Yb) / dTau

Vb = Sqr(VXb ^ 2 + VYb ^ 2)

Xc = FXc(Xb, Yb, Lbc)

Xc1 = FXc(Xb1, Yb1, Lbc)

VXc = (Xc1 - Xc) / dTau

Vc = VXc

Xs = FXs(Xc, Xb, Lcs, Lbs)

Xs1 = FXs(Xc1, Xb1, Lcs, Lbs)

VXs = (Xs1 - Xs) / dTau

Ys = FYs(Yb, Lcs, Lbs)

Ys1 = FYs(Yb1, Lcs, Lbs)

VYs = (Ys1 - Ys) / dTau

Vs = Sqr(VXs ^ 2 + VYs ^ 2)

Psi = Atn(Ys / (Xc - Xs))

Range("F" + LTrim$(Str$(i))) = Xb

Range("G" + LTrim$(Str$(i))) = Yb

Range("H" + LTrim$(Str$(i))) = Vb

Range("I" + LTrim$(Str$(i))) = Xc

Range("J" + LTrim$(Str$(i))) = Vc

Range("K" + LTrim$(Str$(i))) = Xs

Range("L" + LTrim$(Str$(i))) = Ys

Range("M" + LTrim$(Str$(i))) = Vs

i = i + 1

Next Fi1

End Sub

Function FXb(Lab, Omega1, Tau)

FXb = Lab * Cos(Omega1 * Tau)

End Function

Function FYb(Lab, Omega1, Tau)

FYb = Lab * Sin(Omega1 * Tau)

End Function

Function FXc(Xb, Yb, Lbc)

FXc = Xb * Sqr(Lbc ^ 2 - Yb ^ 2)

End Function

Function FXs(Xc, Xb, Lcs, Lbs)

FXs = Xc - Lcs / Lbs * (Xc - Xb)

End Function

Function FYs(Yb, Lcs, Lbs)

FYs = Yb * Lcs / Lbs

End Function

Лабораторная работа 11. Расчет теплообменника.

В приводимой ниже задаче требуется определить распределения температур по длине трубы теплообменника охлаждающей и охлаждаемой жидкостей и стенки трубы. Тип теплообменника «труба в трубе», направление потоков охлаждаемой и охлаждающей жидкостей – навстречу друг другу (противоток) (см. рис.11.1). Внешняя труба теплоизолирована, теплоемкость и теплопроводность материала трубы не учитываем. Не учитываем также передачу тепла вдоль теплонесущих жидкостей.

T - температура, d - диаметр трубы,  - коэффициент теплоотдачи к стенке трубы,  - плотность жидкости, W - скорость течения, Cp - удельная теплоемкость жидкости, L - длина трубы.

Рис.11.1

Для численного решения задачи разобьем трубу на элементарные участки по длине, в которых температуру жидкости будем считать постоянной. Длину элементарного участка примем равной dL. Элементарным объемом жидкости во внутренней трубе будет цилиндр диаметром основания d1 и высотой dL. Элементарным объемом жидкости в наружной трубе будет кольцо наружным диаметром d2, внутренним диаметром d1 и высотой dL.

Тепловой баланс этих элементарных объемов представляют следующие уравнения:

- количество тепла в элементарном объеме внутренней трубы.

- количество тепла, отданного элементарным объемом внутренней трубы в участок стенки трубы длиной dL и температурой Т_S за время контакта с участком стенки (dL/w₁).

- количество тепла в элементарном объеме внутренней трубы.

- количество тепла, полученного элементарным объемом наружной трубы от участка стенки трубы длиной dL и температурой Т_S за время контакта с участком стенки (dL/w₂).

Так как остальными видами передачи тепла пренебрегаем, то

,

откуда следует выражение для определения температуры стенки:

.

Таким образом, зная температуры элементарных объемов жидкостей, можно рассчитать:

1. Температуру стенки Т_S ;

2. Передаваемое количество тепла Q_1охл и Q_2нагр ;

3. Количество тепла элементарных объемов для следующего участка трубы:

Q₁=Q₁-Q_1охл и Q₂=Q₂+ Q_2нагр;

4. Температуры элементарных объемов для следующего участка Т₁ и Т₂.

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

Для противотока известно конечное значение Т₂ (на правом конце трубы), то есть для начала расчета следует задаться каким-либо значением Т₂ для левого конца, а для правого конца трубы расчетное Т₂ должно совпасть с заданным. Для достижения такого совпадения следует организовать итерационный процесс, позволяющий последовательно приблизиться к заданному Т₂. Одним из вариантов такого процесса может быть задание Т₂¹ на левом конце трубы как полусуммы заданных Т₁ и Т₂. В конце расчета, если расчетное Т₂ больше заданного, то в качестве следующего приближения берется полусумма Т₂¹ и заданного Т₂, если расчетное Т₂ меньше заданного, то в качестве следующего приближения берется полусумма Т₂¹ и заданного Т₁. Критерием окончания вычисления берется достижение какой-либо малой разности между расчетным и заданным Т₂.

Пример программы расчета противоточного теплообменника типа «труба в трубе». Итерационный процесс организован выбором конечной температуры охлаждающего потока при заданных начальных температурах охлаждаемого и охлаждающего потоков.

Sub Теплообменник()

T10 = Range("D10")

d1 = Range("D11")

alfa1 = Range("D12")

ro1 = Range("D13")

w1 = Range("D14")

Cp1 = Range("D15")

T20 = Range("E10")

d2 = Range("E11")

alfa2 = Range("E12")

ro2 = Range("E13")

w2 = Range("E14")

Cp2 = Range("E15")

L = Range("D16")

dL = Range("D17")

eps = Range("D18")

m1 = 3.14 * d1 ^ 2 * dL * ro1 / 4

m2 = 3.14 * (d2 ^ 2 - d1 ^ 2) * dL * ro2 / 4

SdL = 3.14 * d1 * dL ^ 2

aw1 = alfa1 / w1

aw2 = alfa2 / w2

T2 = T10

T11 = T10

T21 = T20

j = 1

While Abs(T2 - T20) > eps

T1 = T10

T2 = (T11 + T21) / 2

For i = 1 To L / dL

Ts = (aw1 * T1 + aw2 * T2) / (aw1 + aw2)

Range("F" + LTrim$(Str$(i))) = T1

Range("G" + LTrim$(Str$(i))) = T2

Range("H" + LTrim$(Str$(i))) = Ts

Q1 = Cp1 * m1 * T1

Q1cul = aw1 * (T1 - Ts) * SdL

T1 = (Q1 - Q1cul) / (Cp1 * m1)

Q2 = Cp2 * m2 * T2

Q2hea = aw2 * Abs(Ts - T2) * SdL

T2 = (Q2 - Q2hea) / (Cp2 * m2)

Next i

If T2 < T20 Then

T21 = (T11 + T21) / 2

Else

T11 = (T11 + T21) / 2

End If

Wend

End Sub

Литература

1. Ахо А., Хопкрофт Дж., Ульман Дж.. Построение и анализ вычислительных алгоритмов.- М.:Мир, 1979.

2. Бахвалов Н.С. Численные методы (анализ, алгебра, обыкновенные дифференциальные уравнения). М.: Наука, 1973.

3. Биллиг В.А., Дехтярь М.И., "VBA и Office 97. Офисное программирование" — М.: Издательский отдел "Русская Редакция" ТОО "Channel Trading Ltd." 1998. - 720 с.: ил., компакт-диск.

4. Гарнаев А. Самоучитель VBA. Серия «Самоучитель». BHV-Санкт-Петербург. 1999. – 512с.

5. Кнут Д.. Искусство программирования для ЭВМ. т.3. Сортировка и поиск.- М.: Мир, 1978.

6. Корн Г., Корн Т. Справочник по математике для научных работкников и инженеров. – М.: Наука, 1984. – 831с.

7. Н.Вирт Н. Алгоритмы и структуры данных.- М.:Мир, 1989.

8. Орвис Вильям. Excel для ученых, инженеров и студентов: Пер. с англ. – К.: Юниор, 1999. – 528 с.

9. Островский А.М. Решение уравнений и систем уравнений. М.: ИЛ, 1963.

10. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. – 430с.

11. Сибуя М., Ямамото Т. Алгоритмы обработки данных. - М.:Мир, 1986.

12. Численные методы. Учебник для техникумов. М., «Высш.школа», 1976. – 368с.

Содержание

	Стр.
Введение………………………………………………………………….	3
1. Основы программирования на VBA…………………………………	5
1.1. Создание, редактирование и запись программ……………………	6
1.2. Переменные, константы и типы данных…………………………..	13
1.3. Управляющие конструкции…………………………………………	18
Лабораторная работа 1. Составление алгоритма вычисления значений таблично заданной функции………………………………………	25
Лабораторная работа 2. Поиск экстремума функции…………………..	28
2. Сортировка данных……………………………………………………	31
Лабораторная работа 3. Составление алгоритмов сортировка данных………………………………………………………………………………	36
3. Работа VBA с объектами Excel……………………………………….	39
Лабораторная работа 4. Использование элементов объектного программирования………………………………………………………………	56
4. Численные методы математики……………………………………….	58
4.1. Методы решения нелинейных уравнений…………………………	58
4.2. Численное решение систем линейных алгебраических уравнений	63
4.3. Обработка экспериментальных данных……………………………	72
4.4. Методы численного интегрирования………………………………	81
4.5. Методы численного интегрирования обыкновенных дифференциальных уравнений…………………………………………..	85
4.6. Методы решения линейной краевой задачи для обыкновенных дифференциальных уравнений…………………………………………..	92
Лабораторная работа 5. Определение интеграла функции…………….	96
Лабораторная работа 6. Решение нелинейных уравнений……………....	99
Лабораторная работа 7. Решение систем линейных уравнений……….	101
Лабораторная работа 8. Решение обыкновенного дифференциального уравнения методом Рунге-Кутта………………….	105
Лабораторная работа 9. Построение эмпирической формулы методом наименьших квадратов…………………………………………….	108
Лабораторная работа 10. Моделирование кривошипно-ползунного механизма………………………………………………………………………….	109
Лабораторная работа 11. Расчет теплообменника……………………..	112
Литература	116