Решение уравнений в SciLab

.

Для производной  запишем следующие формулы:

,

,

.

Можем получить три вида разностных уравнений:

,       (2.18)

,       (2.19)

,       (2.20)

.

Разностные уравнения (2.18) аппроксимируют уравнение (2.16) с погрешностью , уравнение (2.19) - с такой же погрешностью, а уравнение (1.7) уже аппроксимирует уравнение (2.16) с погрешностью .

В разностной схеме (2.17) задействованы 4 узла. Конфигурация схемы (2.17) имеет вид:

В схеме (2.19) также участвуют 4 узла, и эта схема имеет вид:

В схеме (2.20) участвуют 5 узлов, и эта схема имеет вид:

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

Для узлов начального (нулевого) слоя значения решения выписываются с помощью начального условия:

Для граничных узлов, лежащих  на прямых и , заменив производные по формулам численного дифференцирования, получаем из граничных условий следующие уравнения:

     (2.21)

Уравнения (2.21) аппроксимируют граничные условия с погрешностью , так как используем односторонние формулы численного дифференцирования. Погрешность аппроксимации можно понизить, если использовать более точные односторонние (с тремя узлами) формулы численного дифференцирования.

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

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

С точки зрения точечной аппроксимации третья схема самая  точная.

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

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

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

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

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

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

 

 

 

 

 

Задача и листингн  программы

Найти распределение температуры  в стержне длиной L, начальная  температура

стержня задается произвольной функцией f(x). Температуры концов стержня равны

u0,t=U1=const , uL,t=U2=const. На боковой поверхности стержня

происходит  теплообмен по закону Ньютона со средой температура которой равна ,

u0 .

//Начальное условие

function y=fi(x)

y=exp(0.15*x)

endfunction

function [u,x,t,r,k]=neiav(N,K,L,T,a,h,U1,U2,u0,eps)

//Функция решения параболического  уравнения методом сеток с

// помощью неявной разностной  схемы. 

//N - количество участков, на которые разбивается интервал по13

х

// (0,L); K - количество участков, на которые разбивается

// интервал по t (0,T);

//a, h - параметры диффренциального уравнения

теплопроводности;

// eps - точность решения СЛАУ) методом Зейделя;

// U1 - температура на левом конце стержня;

// U2 - температура на правом конце стержня;

// Функция neiav возвращает:

// u - матрицу решений в узлах сетки, массив x, массив t,

// r - точность решения системы методом Зейделя,

// к - количество итераций при решении системы методом

Зейделя.

// Вычисляем шаг по х

hx=L/N;

// Вычисляем шаг по  t

delta=T/K;

// Формируем массив x и первый столбец матрицы решений U

// из начального условия

for i=1:N+1

x(i)=(i-1)*hx;

u(i,1)=fi(x(i));

end

// Формируем массив t, первую и последнюю строку матрицы

// решений U из граничных условий

for j=1:K+1

t(j)=(j-1)*delta;

u(1,j)=U1;

u(N+1,j)=U2;

end

// Определяем матрицу  ошибок R и заполняем ее нулями

R(N+1,K+1)=0;

// Вычисляем коэффициент  gamma

gam=a^2*delta/hx^2;

r=1;

k=0;

//Цикл while для организации итерационного процесса при

решении

//системы уравнений методом  Зейделя с точностью eps

while r>eps

// Вычисление матрицы  ошибок R во внутренних точках

//и пересчет значений  u во внутренних точках при решениии СЛАУ

// (12.17) методом Зейделя

for j=2:K+1

for i=2:N

V=gam*(u(i-1,j)+u(i+1,j))/(1+2*gam+delta*hx)+u(i,j-1)/...

(1+2*gam+delta*hx)+delta*h*u0/(1+2*gam+delta*hx);

R(i,j)=abs(V-u(i,j));

u(i,j)=V;

end14

end

// Поиск максимума в  матрице ошибок

r=R(1,1);

for i=1:N+1

for j=1:K+1

if R(i,j)>r

r=R(i,j);

end

end

end

// Увеличение количества  итераций

k=k+1;

end

endfunction

//Вызов функции решения  задачи 

[U,X,T,R,K]=neiav(50,200,5,3,0.4,0.5,1,2.117,30,0.001);

//Построение графика функции

surf(X,T,U');

title('Example 12.2');

xlabel('X');

ylabel('T');

 

 

 

 

 

 

График решения задачи

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Заключение

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

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

 

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

Также изучена  литература, описывающая физическую интерпретацию данного уравнения.

С помощью  среды SciLab были проверены вычисления, проводимые в прикладной программе. Эта проверка показала достаточно высокую точность вычислений прикладной программы. Значения, полученные в среде SciLab, отличаются от вычислений прикладной программы примерно на одну тысячную. Это позволяет судить о достаточно высокой точности вычислений прикладной программы и соответственно о точности вычислений интеграла методом трапеций, который использовался в прикладной программе.