Пусть задана система точек , в которых известны значения функции .
Рассмотрим пример построения интерполяционного многочлена Лагранжа по заданной системе точек (в общем случае для неравноудаленных аргументов). Построим некоторый многочлен таким образом, чтобы его значения совпали со значениями функции, заданными в таблице, для тех же аргументов, то есть . Лагранж предложил строить многочлен степени n в виде:
Здесь в каждом слагаемом коэффициенту ai соответствует разность . Найдем неизвестные коэффициенты , называемые коэффициентами Лагранжа, используя условие .
При :
;
.
Следовательно, коэффициент a1 вычисляется по следующей формуле:
.
При :
;
.
Следовательно, коэффициент a1 вычисляется по формуле:
.
Значения остальных коэффициентов вычисляются аналогично.
С учетом найденных коэффициентов интерполяционный многочлен Лагранжа запишется в виде:
Погрешность формулы (2.16) определяется остаточным членом:
|
,
|
(2.17)
|
где
|
ξ – точка наименьшего промежутка, содержащего все узлы и точку x.
|
Используя условные обозначения суммы и произведения, формулу (2.16) можно записать в сокращенном виде:
В отличие от кусочно-линейной, кусочно-квадратичной интерполяции и кубического сплайна, формула многочлена Лагранжа будет общая для всего интервала , поэтому определять номер интервала для каждой расчетной точки не нужно.
2.6 Кусочно-квадратичная интерполяция функции по таблице значений
Задание: с помощью кусочно-квадратичной интерполяции построить график приближенной функции, проходящей через экспериментальные точки, представленные в таблице 2.1, и вычислить значения неизвестной функции в промежуточных точках.
Таблица 2.1 – Экспериментальные данные для интерполяции
С помощью среды Matlab напишем программу для интерполяции.
Проект будет включать следующие файлы:
Lab4_1.m – скрипт выписывает массивы из таблцы, сортирует их, вычисляет и строит график (листинг 2.1);
Lab4_2.m – скрипт выписывает массивы из таблцы, сортирует их, вычисляет и строит график (листинг 2.2);
Gauss.m – метод Гаусса (листинг 2.3);
Interval.m – интервал (листинг 2.4);
Sort_array.m – сортировка массива (листинг 2.5);
Square_val.m – вычисления (листинг 2.6);
Exchange.m-создает нули функций под диагональю(листинг 2.7);
clc
x=[-8.5 -5 -3 -1.5 -0.5 0 1 2.5 4 5];
y=[8 4 -3 -5.5 -3 -1 1.5 0.5 3 6.5];
x1=[-6 -4 -2.5 -1 1.5 2 3.5 5.5];
[x,y]=sort_point(x,y);
x1=sort(x1);
itr=interval(x,x1);
y1=square_val(x,y,x1,itr);
disp('Значения в расчетных точках');
disp(x1);
disp(y1);
N=10;
n=size(x,2);
x2=zeros(1,1);
x2(1)=x(1);
for i=1:(n-1)
h=(x(i+1)-x(i))/N;
x2i=x(i):h:x(i+1);
x2i=x2i (2:end);
x2=[x2 x2i];
end
itr = interval(x, x2);
y2 = square_val(x, y, x2, itr);
figure;
plot(x, y, 'r*', x1, y1, 'go', x2, y2);
title ('Кусочно - квадратичная интерполяция');
legend ('Эксперементальные точки', 'Точки вычисления', 'График');
xlabel ('x');
ylabel ('y');
grid on;
Листинг 2.1, лист 1 – Главный скрипт
clc
x=[1.5 5.5 8 9 9.75 10.5 13.5 16.5 19];
y=[2 4 7.5 11.5 14 7 3 3 4.5];
x1=[0 4 8.5 9.5 10 12 15];
[x,y]=sort_point(x,y);
x1=sort(x1);
itr=interval(x,x1);
y1=square_val(x,y,x1,itr);
disp('Значения в расчетных точках');
disp(x1);
disp(y1);
N=10;
n=size(x,2);
x2=zeros(1,1);
x2(1)=x(1);
for i=1:(n-1)
h=(x(i+1)-x(i))/N;
x2i=x(i):h:x(i+1);
x2i=x2i (2:end);
x2=[x2 x2i];
end
itr = interval(x, x2);
y2 = square_val(x, y, x2, itr);
figure;
plot(x, y, 'r*', x1, y1, 'go', x2, y2);
title ('Кусочно - квадратичная интерполяция');
legend ('Эксперементальные точки', 'Точки вычисления', 'График');
xlabel ('x');
ylabel ('y');
grid on;
Листинг 2.2, лист 1 – Главный скрипт
function [x, d] = Gauss(A,B)
n = size(A, 1);
C = [A B];
for i = 1:n
[~, index] = max(abs(C(i:n, i)));
index = index + i - 1;
if (index ~= i)
temp = C(i, :);
C(i, :) = C(index, :);
C(index, :) = temp;
end
end
n = size(A);
n = n(1);
for i = 1:n
C(i,1:n + 1) = C(i,1:n + 1) ./ C(i,i);
for j = i + 1:n
C(j,1:n+1) = C(j,1:n+1) - C(i,1:n+1) .* C(j,i);
end;
end;
x = zeros(n,1);
for k = n:-1:1
S = 0;
for i = k + 1:n
S = S + x(i,1) * C(k,i);
end;
x(k,1) = C(k,n + 1) - S;
end;
d = A * x - B;
return
Листинг 2.3, лист 1 – Метод Гаусса
function [itr] = interval(x, x1)
n = length(x);
n1 = length(x1);
itr = zeros(n1, 1);
for i = 1:n1
if (x1(i) < x(1))
itr(i) = 0;
continue;
end
if(x1(i) > x(n))
itr(i) = n;
end
j = 1;
while (j <= n - 1)
if (x1(i) >= x(j) && x1(i) <= x (j+1))
itr(i) = j;
break;
else
j = j+1;
end
end
end
return
end
Листинг 2.4, лист 1 – Интервал
j = 1;
while (j <= n - 1)
if (x1(i) >= x(j) && x1(i) <= x (j+1))
itr(i) = j;
break;
else
j = j+1;
end
end
end
return
end
Листинг 2.5, лист 2 – Интервал
function [x, y] = sort_point(x, y)
n=length(x);
for i=1:(n-1)
for j=1: (n-i)
if (x(j)>x(j+1))
temp=x(j);
x(j)=x(j+1);
x(j+1)=temp;
temp=y(j);
y(j)=y(j+1);
y(j+1)=temp;
end
end
end
end
Листинг 2.6, лист 1 – Сортировка массива
function [y1] = square_val(x, y, x1, itr)
n = length(x);
n1 = length(x1);
a = zeros(n1);
b = zeros(n1);
c = zeros(n1);
y1 = zeros(1, n1);
i = 1;
while (i <= n - 2)
X = [x(i)^2 x(i) 1;
x(i+1)^2 x(i+1) 1;
x(i+2)^2 x(i+2) 1];
Y = [y(i); y(i+1); y(i+2)];
C = Gauss(X, Y);
a(i) = C(1);
b(i) = C(2);
c(i) = C(3);
i = i+1;
end
a(n-1) = a(n-2);
b(n-1) = b(n-2);
c(n-1) = c(n-2);
i = 1;
while (i <= n1)
j = itr(i);
if ((j == 0) || (j == n))
y1(i) = NaN;
end
if ((j > 0) && (j < n))
y1(i) = a(j) * x1(i).^2 + b(j) * x1(i) + c(j);
end
i = i + 1;
end
end
Листинг 2.6, лист 2 – Скрипт вычислений
function z=Exchange(C,i)
k=i+1;
while C(k,i)==0
k=k+1;
end
for j=1:size(C,1)
s=C(i,j);
C(i,j)=C(k,j);
C(k,j)=s;
end
z=C;
end
Листинг 2.7, лист 1 – Функция Exchange
Значения функции в расчетных точках приведены в таблице 2.2.
Таблица 2.2 – Значения функции при кусочно-квадратичной интерполяции
|
Точки вычисления f(x)
|
1
|
x
|
-6
|
-4
|
-2.5
|
-1
|
1.5
|
2
|
3.5
|
5.5
|
y
|
6.2143
|
-0.0238
|
-4.6667
|
-4.5000
|
0.7778
|
0.4444
|
1.8000
|
NaN
|
2
|
x
|
0
|
4
|
8.5
|
9.5
|
10
|
12
|
15
|
|
y
|
NaN
|
2.7308
|
9.5952
|
14.2222
|
11.4000
|
4.5000
|
2.7545
|
|
Графики для первой и второй функции, построенные с помощью программы Matlab, представлены на рисунке 2.3.
а б
Рисунок 2.3 – Квадратичная интерполяция для первой (а) и второй (б) функции
2.7 Интерполяция функции кубическим сплайном по таблице значений
Задание: с помощью интерполяции кубическим сплайном построить график приближенной функции, проходящей через экспериментальные точки, представленные в таблице 2.1, и вычислить значения неизвестной функции в промежуточных точках.
С помощью среды Matlab напишем программу для интерполяции.
Проект будет включать следующие файлы:
spline_main.m – главный скрипт записывающий значения из варианта, вычисляет и строит график (листинг 2.7);
interval.m – интервал (листинг 2.8);
progon.m – прогон (листинг 2.9);
sort_point.m – сортировка точек (листинг 2.10);
spline_val.m – вычисления (листинг 2.11);
variant.m – вариант (листинг 2.12);
[x,y]=sort_point(x,y);
x1=sort(x1);
M=progon(x,y);
itr=interval(x,x1);
y1=spline_val(x,y,x1,itr,M);
disp('Значения в расчетных точках:');
disp(x1);
disp(y1);
a=min([x x1]);
b=max([x x1]);
x2=a:0.1:b;
itr=interval(x,x2);
y2=spline_val(x,y,x2,itr,M);
plot(x,y,'bo',x1,y1,'r*',x2, y2);
title ('Интерполяция кубическим сплайном');
Листинг 2.7, лист 1 – Главный скрипт
legend ('Заданные узлы', 'Расчетные узлы', 'Кривая интерполяции');
xlabel ('x');
ylabel ('y');
grid on;
Листинг 2.7, лист 2 – Главный скрипт
function [itr] = interval(x, x1)
n = length(x);
n1 = length(x1);
itr = zeros(n1, 1);
for i = 1:n1
if (x1(i) < x(1))
itr(i) = 0;
continue;
end
if(x1(i) > x(n))
itr(i) = n;
end
j = 1;
while (j <= n - 1)
if (x1(i) >= x(j) && x1(i) <= x (j+1))
itr(i) = j;
break;
else
j = j+1;
end
end
end
return
Листинг 2.8, лист 1 – Интервал
function [ M ]= progon(x, y)
% входные аргументы: (x,y) - таблица заданных узлов
% выходные аргументы: M - массив вторых производных в узловых точках
n = length(x)-1; % число интервалов
h(1:n)=x(2:n+1)-x(1:n); % длины интервалов
% формирование трехдиагональной СЛАУ
A(1)=0; A(2:n-1)=h(2:n-1); % нижняя диагональ СЛАУ
B(1:n-1)=2*(h(1:n-1)+h(2:n)); % главная диагональ СЛАУ
C(1:n-2)=h(2:n-1); C(n-1)=0; % верхняя диагональ СЛАУ
D=zeros(1:n-1); % правая часть СЛАУ
Листинг 2.9, лист 1 – Прогон
for i=1:n-1
D(i)=6*((y(i+2)-y(i+1))/h(i+1)-(y(i+1)-y(i))/h(i));
end
% прямой ход метода прогонки
Q=zeros(1,n);
R=zeros(1,n);
for i=1:n-1
Q(i+1)=-(C(i)/(B(i)+A(i)*Q(i)));
R(i+1)=(D(i)-A(i)*R(i))/(B(i)+A(i)*Q(i));
end
% обратный ход метода прогонки
M=zeros(1,n-1);
M(n-1)=R(n);
for i=n-2:-1:1
M(i)=Q(i+1)*M(i+1)+R(i+1);
end
% дополняем массив произодных краевыми условиями M(a)=M(b)=0
M = [0, M, 0];
End
Листинг 2.9, лист 2 – Прогон
function [x,y] = sort_point(x,y)
n=length(x);
for i=1:(n-1)
flag=0;
for j=1: (n-i)
if (x(j)>x(j+1))
temp=x(j); x(j)=x(j+1); x(j+1)=temp;
temp=y(j); y(j)=y(j+1); y(j+1)=temp;
flag=1;
end
end
if flag==0
break;
end
end
Листинг 2.10, лист 1 – Сортировка точек
function [y1]=spline_val(x,y,x1,itr,M)
n=length(x)-1;
n1=length(x1);
y1=zeros(1,n1);
h(1:n)=x(2:n+1)-x(1:n);
i=1;
while(i<=n1)
j=itr(i);
if(j==0)
y1(i)=y(1)+((x(1)-x(2))*M(2)/6+(y(2)-y(1))/(x(2)-x(1)))*(x1(i)-x(1));
i=i+1;
end
if(j>0 && j<=n)
y1(i)=(1/(6*h(j)))*((M(j)*(x(j+1)- x1(i))^3)+M(j+1)*(x1(i)-x(j))^3)+(1/h(j))*((y(j)-((M(j)*h(j)^2)/6))*(x(j+1)-x1(i))+(y(j+1)-((M(j+1)*h(j)^2)/6))*(x1(i)-(x(j))));
i=i+1;
end
if(j==n+1)
y1(i)=y(n+1)+((x(n+1)-x(n))*M(n)/6+(y(n+1)-y(n))/(x(n+1)-x(n)))*(x1(i)-x(n+1));
i=i+1;
Листинг 2.11, лист 1 – Вычисления
end
end
end
Листинг 2.11, лист 2 – Вычисления
function [ x,y,x1 ] = variant( var,table )
% Наборы экспериментальных и расчетных точек для интерполяции
% var - номер варианта
% table - номер таблицы (1 или 2)
switch(var)
case 2
if table==1
x=[-9 -7.5 -5.5 -3.5 -0.5 1.5 3.5 4.5 7 9];
y=[8 2.5 -2.5 -3.75 -1 3.5 6.5 4.5 6 6];
x1=[-8.5 -6 -2 0 1 2.5 5.5 8];
elseif table==2
x=[1.5 5.5 8 9 9.75 10.5 13.5 16.5 19];
y=[2 4 7.5 11.5 14 7 3 3 4.5];
x1=[0 4 8.5 9.5 10 12 15];
end
end
end
Листинг 2.12, лист 1 – Наш вариант
Значения функции в расчетных точках приведены в таблице 2.3.
Таблица 2.3 – Значения функции при интерполяции кубическим сплайном
|
Точки вычисления f(x)
|
1
|
x
|
-6
|
-4
|
-2.5
|
-1
|
1.5
|
2
|
3.5
|
5.5
|
y
|
6.2005
|
0.6021
|
-4.4604
|
-4.6362
|
1.4271
|
0.8864
|
1.6454
|
8.3557
|
2
|
x
|
0
|
4
|
8.5
|
9.5
|
10
|
12
|
15
|
|
y
|
1.3728
|
3.1253
|
9.1147
|
14.0456
|
12.3930
|
0.4735
|
3.6721
|
|
Графики для первой и второй функции, построенные с помощью программы Matlab, представлены на рисунке 2.4.
а б
Рисунок 2.4 – Интерполяция сплайном для первой (а) и второй (б) функции
2.8 Интерполяция функции многочленом Лагранжа по таблице значений
Задание: с помощью интерполяции многочленом Лагранжа построить график приближенной функции, проходящей через экспериментальные точки, представленные в таблице 2.1, и вычислить значения неизвестной функции в промежуточных точках.
С помощью среды Matlab напишем программу для интерполяции.
Проект будет включать следующие файлы:
lagr_main.m – главный скрипт выписывающий значения из варианта, вычисляющий и строящий график (листинг 2.13);
lagr_val.m – вычисления лагранжа (листинг 2.14);
sort_point.m – сортировка точек (листинг 2.15);
variant.m - вариант (листинг 2.16);
Значения функции в расчетных точках приведены в таблице 2.4.
[x,y]=sort_point(x,y);
x1=sort(x1);
y1=lagr_val(x,y,x1);
disp('Значения в расчетных точках:');
disp(x1);
disp(y1);
a=min(x);
b=max(x);
x2=a:0.1:b;
y2=lagr_val(x,y,x2);
plot(x,y,'bo',x1,y1,'r*',x2, y2);
title ('Многочлен Лагранжа');
legend ('Заданные узлы', 'Расчетные узлы', 'Кривая интерполяции');
xlabel ('x');
ylabel ('y');
grid on;
Листинг 2.13, лист 1 – Главный скрипт
function [y1]=lagr_val(x,y,x1)
n=length(x);
n1=length(x1);
y1=zeros(1,n1);
k=1;
while(k<=n1)
if((x1(k)x(n)))
y1(k)=0;
k=k+1;
else
i=1;
while(i<=n)
P=1;
j=1;
while(j<=n)
if(i~=j)
P=P*(x1(k)-x(j))/(x(i)-x(j));
end
j=j+1;
end
y1(k)=y1(k)+P*y(i);
i=i+1;
end
k=k+1;
end
end
end
Листинг 2.14, лист 1 – Вычисления лагранжа
function [x,y] = sort_point(x,y)
n=length(x);
for i=1:(n-1)
flag=0;
for j=1: (n-i)
if (x(j)>x(j+1))
temp=x(j); x(j)=x(j+1); x(j+1)=temp;
temp=y(j); y(j)=y(j+1); y(j+1)=temp;
flag=1;
end
end
if flag==0
break;
end
end
Листинг 2.15, лист 1 – Сортировка точек
function [ x,y,x1 ] = variant( var,table )
% Наборы экспериментальных и расчетных точек для интерполяции
% var - номер варианта
% table - номер таблицы (1 или 2)
switch(var)
case 13
if table==1
x=[-8.5 -5 -3 -1.5 -0.5 0 1 2.5 4 5];
y=[8 4 -3 -5.5 -3 -1 1.5 0.5 3 6.5];
x1=[-6 -4 -2.5 -1 1.5 2 3.5 5.5];
elseif table==2
x=[1.5 5.5 8 9 9.75 10.5 13.5 16.5 19];
y=[2 4 7.5 11.5 14 7 3 3 4.5];
x1=[0 4 8.5 9.5 10 12 15];
end
end
Листинг 2.16, лист 1 – Вариант
Таблица 2.4 – Значения функции при интерполяции многочленом Лагранжа
|
Точки вычисления f(x)
|
1
|
x
|
-6
|
-4
|
-2.5
|
-1
|
1.5
|
2
|
3.5
|
5.5
|
y
|
48.2638
|
-2.1426
|
-4.1887
|
-4.6868
|
1.4688
|
0.9274
|
1.7147
|
0
|
2
|
x
|
0
|
4
|
8.5
|
9.5
|
10
|
12
|
15
|
|
y
|
0
|
-318.59
|
8.3240
|
13.9518
|
12.9194
|
-24.8923
|
104.56
|
|
Графики для первой и второй функции, построенные с помощью программы Matlab, представлены на рисунке 2.5.
а б
Рисунок 2.5 – Многочлен Лагранжа для первой (а) и второй (б) функции
Вывод
Если сравнивать три вида интерполяции, рассмотренные выше, то можно сделать вывод, что для заданного набора значений кубический сплайн – это самый точный метод интерполяции. Использование многочленов третьей степени и согласование значений производных в узлах позволяет получить гладкую кривую без значительных отклонений.
Квадратичная интерполяция основана на уравнении параболы. Она не учитывает значения производных, поэтому дает более высокую погрешность. Кроме того, график функции получается ломаным.
Интерполяция многочленом Лагранжа является более простой в понимании и реализации, потому что используется одна форума на всем интервале интерполяции. Однако, повышение степени многочлена приводит к нежелательным колебаниям графика, которые особенно заметны на краях.
3 Решение обыкновенных дифференциальных уравнений
Дифференциальные уравнения являются основным математическим инструментом моделирования и анализа разнообразных явлений и процессов в науке и технике.
Методы их решения подразделяются на два класса:
аналитические методы, в которых решение получается в виде аналитических функций;
численные (приближенные) методы, где искомые интегральные кривые получают в виде таблиц их численных значений.
Применение аналитических методов позволяет исследовать полученные решения методами математического анализа и сделать соответствующие выводы о свойствах моделируемого явления или процесса. К сожалению, с помощью таких методов можно решать достаточно ограниченный круг реальных задач. Численные методы позволяют получить с определенной точностью приближенное решение практически любой задачи.
Решить дифференциальное уравнение
численным методом означает, что для заданной последовательности аргументов и числа , не определяя аналитического вида функции , найти значения , удовлетворяющие начальным условиям:
.
3.1 Метод Эйлера
Этот метод является сравнительно грубым и применяется в основном для ориентировочных расчетов. Однако идеи, положенные в основу метода Эйлера, являются исходными для ряда других численных методов.
Пусть дано обыкновенное дифференциальное уравнение (ОДУ) с начальными условиями (задача Коши):
и удовлетворяются условия существования и единственности решения.
Требуется найти решение задачи Коши (3.2) на отрезке . Найдем решение в виде таблицы . Для этого разобьем отрезок на n равных частей и построим последовательность
,
где – шаг интегрирования.
Проинтегрируем исходное уравнение на отрезке :
Полученное соотношение можно переписать как
Если считать подынтегральную функцию постоянной на участке и равной значению в начальной точке этого интервала , то получим
.
Подставляя полученный результат в формулу (3.3), получим основную расчетную формулу метода Эйлера:
Так как для вычисления каждого последующего значения yk+1 требуется только одно значение yk, полученное на предыдущем шаге, то метод Эйлера относят к одношаговым методам.
Вычисление значений осуществляется с использованием формулы (3.4) следующим образом. По заданным начальным условиям и y0, полагая в выражении (3.4), вычисляется значение
Далее, определяя значение аргумента x по формуле , используя найденное значение y1 и полагая в формуле (3.4) , вычисляем следующее приближенное значение интегральной кривой как
Поступая аналогичным образом при , определяем все остальные значения yk, в том числе последнее значение
которое соответствует значению аргумента .
Погрешность метода Эйлера составляет .
Вычислив приближенные значения неизвестной функции в точках , для более гладкого представления искомой интегральной кривой можно воспользоваться одним из методов интерполяции, например, кубическим сплайном.
3.2 Метод Эйлера-Коши
Метод Эйлера-Коши, как и классический метод Эйлера, относится к одношаговым методам, но каждый шаг включает два этапа.
На первом этапе выполняется приблизительная оценка значения функции по обычной формуле Эйлера:
На втором этапе полученное значение используется для более точного расчета по формуле:
Погрешность метода Эйлера составляет , что на порядок выше, чем у метода Эйлера.
3.3 Решение ОДУ методом Эйлера
Задание: найдите методом Эйлера численное решение ОДУ первого порядка на отрезке от a = 0 до b = 2 с шагом интегрирования h = 0,2:
Получите таблицу приближенных значений неизвестной функции y(x) на отрезке интегрирования с постоянным шагом h.
Аппроксимируйте полученные табличные значения кубическим сплайном и постройте приближенный график функции y(x) на отрезке интегрирования, используя скрипты из главы 2. При построении графика шаг интерполяции примите равным 0,1·h.
С помощью среды Matlab напишем программу для интерполяции.
Проект будет включать следующие файлы:
Euler.m – основной скрипт в который записываются исходные данные (листинг 3.1);
euler_main.m – главный скрипт с построением графика и вывода значений(листинг 3.2);
interval.m – скрипт задающий интервал(листинг 3.3);
function [y]=Euler(x,h)
n=length(x);
y=zeros(1, n);
y(1)=0.599;
i=1;
Листинг 3.1, лист 1 – Скрипт с исходными данными
for x=0:h:2
if (iy(i+1)=y(i)+h*(3*(1/.1+x^2 - 1/.x+4 -atanx + log(x+4)-2)*y(i));
i=i+1;
end
end
end
Листинг 3.1, лист 2 – Скрипт с исходными данными
h=0.2;
x=-2 :h:0;
y=Euler(x,h);
disp('Значения x');
disp(x);
disp('Значения y');
disp(y);
M=progon(x,y);
a=min(x);
b=max(x);
x2=a:0.1*h:b;
itr=interval(x,x2);
y2=spline_val(x,y,x2,itr,M);
figure(1);
plot(x,y,'*r',x2,y2);
title ('Решение ОДУ методом Эйлера');
legend ('Расчетные точки', 'Интерполяционная функция');
xlabel ('x');
ylabel ('y');
grid on;
Листинг 3.2, лист 1 – Главный скрипт
function [itr] = interval(x, x1)
n = length(x);
n1 = length(x1);
itr = zeros(n1, 1);
for i = 1: n1
if (x1(i) < x(1))
itr(i) = 0;
continue;
end
if (x1(i) > x(n))
itr(i) = n;
end
j = 1;
while (j <= n - 1)
if (x1(i) >= x(j) && x1(i) <= x (j+1))
itr(i) = j;
break;
else
j = j+1;
end
end
end
end
Листинг 3.3, лист 1 – Интервал
Численное решение ОДУ приведено в таблице 3.1.
Таблица 3.1 – Численное решение ОДУ методом Эйлера
K
|
1
|
2
|
3
|
4
|
5
|
6
|
7
|
8
|
9
|
10
|
11
|
|
Достарыңызбен бөлісу: |