Ошибки в разработке программы для оценки методов решения дифференциальных уравнений и построения графиков
11.07.2019, 15:51. Показов 715. Ответов 1
Добрый день! Моя программа состоит из следующих подпрограмм:
- подпрограмма определения функции для первой производной;
- подпрограмма метода Эйлера;
- подпрограмма исправленного метода Эйлера;
- подпрограмма модифицированного метода Эйлера;
- подпрограмма метода Рунге-Кутта 4-го порядка;
- подпрограмма оценки результатов и построения графиков.
В зависимости от команды, поступившей в оператор switch, выбирается метод решения уравнения и его порядок для оценки результата, а также строится график для решения с шагами h и h/2. При оценке результатов я брал в качестве основного интервал с большим шагом и считал разность между узлами с номером i у первого интервала и номером 2i – 1 у второго интервала с половинчатым шагом. Я использовал элементы 2i – 1, так как нумерация в массивах MatLab начинается не с нуля, а с единицы. После этого я выбирал максимальную погрешность и строил графики для шагов h и h/2. При попытке запустить программу, я увидел предупреждение, что индекс превысил размерность матрицы, и ошибку в строке 8. Я читал, что размерность массивов значений для двух кривых у функции plot должна быть одинаковой, в противном случае графики нужно разносить по разным окнам. Я убрал subplot, но ошибки остались. Если в первом случае всё понятно (надо подгонять размерности для аргумента x и функции y), то почему компилятор ругается на массив оценок r из линии 8? Покажите, пожалуйста, где и на что нужно исправить. И вообще, функция plot идёт после той злополучной строки. Тогда где ещё могло возникнуть превышение размера массива?
Функция для правой части уравнения:
| Matlab M | 1
2
| function d_y = func(x,y) % функция первой производной
d_y = (cos(x)/(x + 1)) - 0.5*y^2; |
|
Метод Эйлера:
| Matlab M | 1
2
3
4
5
6
7
| function [y, x] = eiler(x0,b,h,y0) % метод Эйлера
x = x0:h:b; % определение переменной x, принимающей значения от x0 до b с шагом h
y = ones(1,length(x)+1);
y(1) = y0; % определение начального значения переменной y
for i=2:length(y)
y(i) = y(i-1) + h*func(x(i-1),y(i-1)); % определение остальных значений переменной y
end |
|
Исправленный метод Эйлера:
| Matlab M | 1
2
3
4
5
6
7
| function [y, x] = eiler_2(x0,b,h,y0) % исправленный метод Эйлера
x = x0:h:b; % определение переменной x, принимающей значения от x0 до b с шагом h
y = ones(1,length(x)+1);
y(1) = y0; % определение начального значения переменной y
for i=2:length(y)
y(i) = y(i-1) + (h*(func(x(i-1),y(i-1)) + func(x(i),y(i-1) + h*func(x(i-1),y(i-1)))))/2; % определение остальных значений переменной y
end |
|
Модифицированный метод Эйлера:
| Matlab M | 1
2
3
4
5
6
7
| function [y, x] = eiler_modify(x0,b,h,y0) % модифицированный метод Эйлера
x = x0:h:b; % определение переменной x, принимающей значения от x0 до b с шагом h
y = ones(1,length(x)+1);
y(1) = y0; % определение начального значения переменной y
for i=2:length(y)
y(i) = y(i-1) + h*func(x(i-1) + (h/2),y(i-1) + (h/2)*func(x(i-1),y(i-1))); % определение остальных значений переменной y
end |
|
Метод Рунге-Кутта 4-го порядка:
| Matlab M | 1
2
3
4
5
6
7
8
9
10
11
| function [y, x] = runge_kutta_4(x0,b,h,y0) % метод Рунге-Кутта 4-го порядка
x = x0:h:b; % определение переменной x, принимающей значения от x0 до b с шагом h
y = ones(1,length(x)+1);
y(1) = y0; % определение начального значения переменной y
for i=2:length(y)
k1 = h*func(x(i-1),y(i-1));
k2 = h*func(x(i-1) + (h/2),y(i-1) + (k1/2));
k3 = h*func(x(i-1) + (h/2),y(i-1) + (k2/2));
k4 = h*func(x(i-1) + h,y(i-1) + k3);
y(i) = y(i-1) + ((k1 + 2*k2 + 2*k3 + k4)/6); % определение остальных значений переменной y
end |
|
Функция для оценки результатов и построения графиков:
| Matlab M | 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
| function plot_and_estimation(x0,b,h,y0,select_method) % в зависимости от выбранного метода оценить погрешность и построить график
switch select_method
case 'eiler'
[y_1, x_1] = eiler(x0,b,h,y0);
[y_2, x_2] = eiler(x0,b,h/2,y0);
r = ones(1,length(y_1));
for i=1:length(r)
r(i) = abs(y_1(i) - y_2(2*i-1));
end
r_max = r(1);
for k=1:length(r)
if (r_max < r(k))
r_max = r(k);
end
end
disp(['Погрешность метода Эйлера равна ' num2str(r_max)]);
subplot(2,1,1);
plot(x_1,y_1);
title('Метод Эйлера с шагом h');
subplot(2,1,2);
plot(x_2,y_2);
title('Метод Эйлера с шагом h/2');
case 'eiler_correct'
[y_1, x_1] = eiler_2(x0,b,h,y0);
[y_2, x_2] = eiler_2(x0,b,h/2,y0);
r = ones(1,length(y_1));
for i=1:length(r)
r(i) = abs(y_1(i) - y_2(2*i-1))/3;
end
r_max = r(1);
for k=1:length(r)
if (r_max < r(k))
r_max = r(k);
end
end
disp(['Погрешность исправленного метода Эйлера равна ' num2str(r_max)]);
subplot(2,1,1);
plot(x_1,y_1);
title('Исправленный метод Эйлера с шагом h');
subplot(2,1,2);
plot(x_2,y_2);
title('Исправленный метод Эйлера с шагом h/2');
case 'eiler_modify'
[y_1, x_1] = eiler_modify(x0,b,h,y0);
[y_2, x_2] = eiler_modify(x0,b,h/2,y0);
r = ones(1,length(y_1));
for i=1:length(r)
r(i) = abs(y_1(i) - y_2(2*i-1))/3;
end
r_max = r(1);
for k=1:length(r)
if (r_max < r(k))
r_max = r(k);
end
end
disp(['Погрешность модифицированного метода Эйлера равна ' num2str(r_max)]);
subplot(2,1,1);
plot(x_1,y_1);
title('Модифицированный метод Эйлера с шагом h');
subplot(2,1,2);
plot(x_2,y_2);
title('Модифицированный метод Эйлера с шагом h/2');
case 'runge_kutta_4'
[y_1, x_1] = runge_kutta_4(x0,b,h,y0);
[y_2, x_2] = runge_kutta_4(x0,b,h/2,y0);
r = ones(1,length(y_1));
for i=1:length(r)
r(i) = abs(y_1(i) - y_2(1 + 2*(i-1)))/15;
end
r_max = r(1);
for k=1:length(r)
if (r_max < r(k))
r_max = r(k);
end
end
disp(['Погрешность метода Рунге-Кутта 4-го порядка равна ' num2str(r_max)]);
subplot(2,1,1);
plot(x_1,y_1);
title('Метод Рунге-Кутта 4-го порядка с шагом h');
subplot(2,1,2);
plot(x_2,y_2);
title('Метод Рунге-Кутта 4-го порядка с шагом h/2');
otherwise
disp('Неверный ввод команды!')
end |
|
Результат работы программы:
| Matlab M | 1
2
3
4
5
| >> plot_and_estimation(0,1,0.1,0.2,'eiler');
Index exceeds matrix dimensions.
Error in plot_and_estimation (line 8)
r(i) = abs(y_1(i) - y_2(2*i-1)); |
|
В изображении пример нумерации для массивов, начинающихся с единицы.
0
|
(
Сообщение от БузинВладимир
