Решение системы дифференциальных уравнений Методом Рунге-Кутта 4 порядка

@Jack Frost · Регистрация: 26.12.2012

Студворк — интернет-сервис помощи студентам

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

Система дифф. уравнений:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\frac{dM}{d\xi}=\frac{{R}_{M}}{\xi nR}\\ \frac{d\delta }{d\xi }=\frac{\delta {R}_{\delta }}{n\xi (M-\xi )R}\\ \frac{d\Pi }{d\xi }=\frac{\delta (\xi -M){R}_{\Pi }}{n\xi R}\\$

Где:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?R=\gamma \Pi -\delta {(M-\xi )}^{2}\\ {R}_{M}=-2(n-1)\xi \Pi -2\gamma nM\Pi +(n-1)\xi \delta M(M-\xi )\\ {R}_{\delta }={R}_{M}-2nMR\\ {R}_{\Pi }={R}_{M}-\frac{(n-1)\xi MR}{M-\xi }\\ $
Вводные:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\gamma =\frac{5}{3}\\ n=0.688377\\ \xi=1\leq \xi \lt \infty\\$
M = 0.75
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\delta\\$ = 4;
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\Pi\\$ = 0.75

Matlab M

function dy=fun(ksi,y)
dy = zeros(3,1);   % вектор колонка с нулевыми элементами
n=0.688377;
gam=5./3;
R=((gam.*y(3))-y(2).*((y(1)-ksi).^2));
Rm=-2.*(n-1).*ksi.*y(3)-2.*gam.*n.*y(1).*y(3)+(n-1).*ksi.*y(2).*y(1).*(y(1)-ksi);
Rdel=Rm-2.*n.*y(1).*R;
Rp=Rm-((n-1).*ksi.*y(1).*R./(y(1)-ksi));
%Запишем систему ДУ
dy(1) =(Rm./(ksi.*n.*R));
dy(2) =(y(2).*Rdel)./(n.*ksi.*(y(1)-ksi).*R);
dy(3)=(y(2).*(ksi-y(1)).*Rp)./(n.*ksi.*R);

Matlab M

clear all; clc;
ksi=[1 10];
y0(1)=0.75;
y0(2)=4;
y0(3)=0.75;
[T,Y]=ode23(@fun,ksi,y0);
figure
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'--')
grid on
title('График искомых функций');
legend('M','Del','P');
xlabel('ksi');
disp('Вывод значений функций')
[Y(:,1),Y(:,2),Y(:,3)]

Нужно решить систему Методом Рунге-Кутты 4 порядка.
Просмотрел много тем на форуме и нашел несколько похожих:
1. Решение ОДУ методом Рунге-Кутта
2. Метод Рунге-Кутты для решения системы ОДУ
3. Решить систему ДУ методом Рунге-Кутта с апостериорным рассчетом шага

Из этого возникли вопросы:
Как именно задать R, Rm, Rdel, Rp и начальные условия для gam, n, ksi и M, Delta, P?

Добавлено через 1 час 9 минут
Если пробовать привести к виду из Метод Рунге-Кутты для решения системы ОДУ, то

Matlab M

function dy=five(ksi,y,a,b,M)
dy = zeros(3,1);    % вектор колонка с нулевыми элементами
n=0.688377;
gam=5./3;
R1=(n-1).*ksi.*((y(1)-ksi).*y(2).*y(1)-2.*y(3));
R2=n.*ksi.*(gam.*y(3)-((y(1)-ksi).^2).*y(2));
% далее запишем систем ДУ
dy(1) =(R1-2.*gam.*n.*y(1).*y(3))./R2;
dy(2) =(y(2).*(2.*y(1).*y(2).*n.*((y(1)-ksi).^2)-R1))./(R2.*(y(1)-ksi));
dy(3)=((y(2).*y(3))./R2).*(2.*(n.*gam.*y(1)+ksi.*(n-1)).*(y(1)-ksi)-(n-1).*gam.*ksi.*y(1));

Matlab M

clear all; clc;
ksi=[1 10];
y0(1)=0.75;
y0(2)=4;
y0(3)=0.75;
[T,Y]=ode23(@five,ksi,y0);
figure
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'--')
grid on
title('График искомых функций');
legend('M','Del','P');
xlabel('ksi');
disp('Вывод значений функций')
[Y(:,1),Y(:,2),Y(:,3)];

Matlab M

clear all; clc;
% Программа решения методом Рунге-Кутты 4-го порядка
% системы дифференциальных уравнений
% y1' = f1(x,y1,y2)
% y2' = f2(x,y1,y2)
% с начальными условиями y1(0)=y10 и y2(0)=y20
% на интервале [a,b] с перемеенным шагом h
% и точностью eps
% функции правой части
n=0.688377;
gam=5./3;
t=[1 10]; % ksi=t
R1=(n-1).*t.*((x-t).*y.*x-2.*z);
R2=n.*t.*(gam.*z-((x-t).^2).*y);
f1 = @(t,x,y,z) (R1-2.*gam.*n.*x.*z)./R2;
f2 = @(t,x,y,z) (y.*(2.*x.*y.*n.*((x-t).^2)-R1))./(R2.*(x-t));
f3 = @(t,x,y,z) ((y.*z)./R2).*(2.*(n.*gam.*x+t.*(n-1)).*(x-t)-(n-1).*gam.*t.*x);
a = 1; % нижняя граница диапазона
b = 10; % верхняя граница диапазона
eps = 0.01; % точность
h = (b-a)/100; % начальное значение шага
T=a;
X(1)=0.75;
Y(1)=4;
Z(1)=0.75;
% решение системы методом Рунге-Кутты 4-го порядка
i = 1; k = 0;
while T(end) <= b
    % решение первого уравнения системы
    while 1
        % считаем с исходным шагом
        K1 = h*f1( T(i),X(i),Y(i),Z(i));
        L1 = h*f2( T(i),X(i),Y(i),Z(i));
        B1=  h*f3( T(i),X(i),Y(i),Z(i));
        K2 = h*f1( T(i)+h/2, X(i)+K1/2,Y(i)+L1/2, Z(i)+B1/2);
        L2 = h*f2( T(i)+h/2, X(i)+K1/2,Y(i)+L1/2, Z(i)+B1/2);
        B2 = h*f3( T(i)+h/2, X(i)+K1/2,Y(i)+L1/2, Z(i)+B1/2);
        K3 = h*f1( T(i)+h/2, X(i)+K2/2,Y(i)+L2/2, Z(i)+B2/2);
        L3 = h*f2( T(i)+h/2, X(i)+K2/2,Y(i)+L2/2, Z(i)+B2/2);
        B3 = h*f3( T(i)+h/2, X(i)+K2/2,Y(i)+L2/2, Z(i)+B2/2);
        K4 = h*f1( T(i)+h, X(i)+K3, Y(i)+L3, Z(i)+B3);
        L4 = h*f2( T(i)+h, X(i)+K3, Y(i)+L3, Z(i)+B3);
        B4 = h*f3( T(i)+h, X(i)+K3, Y(i)+L3, Z(i)+B3);
        delta11 = X(i) + (K1 + 2*K2 + 2*K3 + K4)/6;
        delta12 = Y(i) + (L1 + 2*L2 + 2*L3 + L4)/6;
        delta13 = Z(i) + (B1 + 2*B2 + 2*B2 + B4)/6;
        
        h = h/2; % уменьшаем шаг в двое
        
        K1 = h*f1( T(i),X(i),Y(i),Z(i));
        L1 = h*f2( T(i),X(i),Y(i),Z(i));
        B1=  h*f3( T(i),X(i),Y(i),Z(i));
        K2 = h*f1( T(i)+h/2, X(i)+K1/2,Y(i)+L1/2, Z(i)+B1/2);
        L2 = h*f2( T(i)+h/2, X(i)+K1/2,Y(i)+L1/2, Z(i)+B1/2);
        B2 = h*f3( T(i)+h/2, X(i)+K1/2,Y(i)+L1/2, Z(i)+B1/2);
        K3 = h*f1( T(i)+h/2, X(i)+K2/2,Y(i)+L2/2, Z(i)+B2/2);
        L3 = h*f2( T(i)+h/2, X(i)+K2/2,Y(i)+L2/2, Z(i)+B2/2);
        B3 = h*f3( T(i)+h/2, X(i)+K2/2,Y(i)+L2/2, Z(i)+B2/2);
        K4 = h*f1( T(i)+h, X(i)+K3, Y(i)+L3, Z(i)+B3);
        L4 = h*f2( T(i)+h, X(i)+K3, Y(i)+L3, Z(i)+B3);
        B4 = h*f3( T(i)+h, X(i)+K3, Y(i)+L3, Z(i)+B3);
        delta21 = X(i) + (K1 + 2*K2 + 2*K3 + K4)/6;
        delta22 = Y(i) + (L1 + 2*L2 + 2*L3 + L4)/6;
        delta23 = Z(i) + (B1 + 2*B2 + 2*B2 + B4)/6;
        
        if max(abs([delta11-delta21, delta12-delta22,delta13-delta23])) > eps % Если разность больше точности
            h = h/2; % уменьшаем шаг
        else % если подходит
            h = 2*h; % возвращаем назад
            % берем первые значения приращения
            T(i+1) = T(i) + h;
            X(i+1) = delta11;
            Y(i+1) = delta12;
            Z(i+1) = delta13;
            break
        end
    end
    i = i+1;
    h = 2*h; % увеличиваем шаг для следующей итерации
    k = k+1; % счетчик итераций
    if k > 50000
        disp('Решение не найдено (превышен лимит итераций)')
        break
    end
end
 
plot(T,X,T,Y,T,Z) % построение графика
grid on % наложение сетки
legend('X(T)','Y(T)','Z(T)'); % подпись кривых
title('Метод Рунге-Кутты 4-го порядка') % заголовок окна
xlabel('X') % надпись по оси X
ylabel('Y') % надпись по оси Y

И получается проблема в

Matlab M

n=0.688377;
gam=5./3;
t=[1 10]; % ksi=t
R1=(n-1).*t.*((x-t).*y.*x-2.*z);
R2=n.*t.*(gam.*z-((x-t).^2).*y);

@SSC · 15.06.2017, 07:02

У Вас есть функция

Matlab M

1	function dy=five(ksi,y,a,b,M)

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

@Jack Frost · 15.06.2017, 10:18 **[ТС]**

Сообщение от SSC

У Вас есть функция
Matlab MВыделить код
1
function dy=five(ksi,y,a,b,M)
которая находит производные в заданной точке.
И нафига Вам сдались функции f1, f2, f3.
Вы повторяете теже проблемы как в предыдущей теме (как говорится лыко-мочало, начинай сначала)
Посмотрите свежую тему, там все сделано для одного уравнения, для трех отличается только тем, что из функции получаете не одно значение, а массив у, и обрабатывайте как массив.
Построить графики решений обыкновенных дифференциальных уравнений

Вы имели в виду это?

Matlab M

%Метод Рунге-Кутта
y4(1)=y1(1); 
b=0.5;
n=(b-x(1))/h+1;
for i=2:n
    r1=h*dx(x(i-1),y4(i-1));
    r2=h*dx((x(i-1)+(h/2)),(y4(i-1)+(r1/2)));
    r3=h*dx((x(i-1)+(h/2)),(y4(i-1)+(r2/2)));
    r4=h*dx((x(i-1)+h),(y4(i-1)+r3));
    y4(i)=y4(i-1)+((r1+(2*r2)+(2*r3)+r4)/6);
    x(i)=x(i-1)+h;
end
subplot(2,2,4);
plot(x,y4);
grid on
title('Метод Рунге-Кутта');

Т.е. обращение будет в виде

Matlab M

1	r1=h*@five(x(i-1),y4(i-1));

?

@SSC · 15.06.2017, 10:40

Сообщение от Jack Frost

обращение будет в виде

Почти.
Т.к. y4 это будет двумерный массив, то обращение вот так

Matlab M

1	r1=h*@five(x(i-1),y4(:,i-1));

и 10 стр.

Matlab M

1	y4(:,i)=y4(:,i-1)+((r1+(2r2)+(2r3)+r4)./6);

и начальные условия надо как-то так задавать

Matlab M

1
2
3

y4(1,1)=0;
y4(2,1)=0;
y4(3,1)=0;

@Jack Frost · 15.06.2017, 10:53 **[ТС]**

Сообщение от SSC

Почти.
Т.к. y4 это будет двумерный массив, то обращение вот так

Matlab M

1	r1=h*@five(x(i-1),y4(:,i-1));

А почему именно?

Matlab M

x(i-1)

Сообщение от SSC

и начальные условия надо как-то так задавать

Matlab M

1
2
3

y4(1,1)=0;
y4(2,1)=0;
y4(3,1)=0;

Спасибо!

Т.е.

Matlab M

n=0.688377;
gam=5./3;
R1=(n-1).*ksi.*((y(1)-ksi).*y(2).*y(1)-2.*y(3));
R2=n.*ksi.*(gam.*y(3)-((y(1)-ksi).^2).*y(2));

из первого файла мне не нужно, т.к. все есть в функции.

А

Matlab M

1
2
3

y0(1)=0.75;
y0(2)=4;
y0(3)=0.75;

Задается через y4(1,1); y4(2,1); y4(3,1); ?

@SSC · 15.06.2017, 11:19

Сообщение от Jack Frost

из первого файла мне не нужно, т.к. все есть в функции.

Да.

Сообщение от Jack Frost

Задается через y4(1,1); y4(2,1); y4(3,1); ?

Да.

Сообщение от Jack Frost

А почему именно?

Определение производной происходит в предыдущей точке расчета.
Обратите внимание на цикл от 2 до n

@Jack Frost · 15.06.2017, 11:28 **[ТС]**

Сообщение от SSC

Определение производной происходит в предыдущей точке расчета.
Обратите внимание на цикл от 2 до n

Спасибо за ответ!
По идее, должно получаться так
five.m

Matlab M

function dy=five(ksi,y)
dy = zeros(3,1);    % вектор колонка с нулевыми элементами
n=0.688377;
gam=5./3;
R1=(n-1).*ksi.*((y(1)-ksi).*y(2).*y(1)-2.*y(3));
R2=n.*ksi.*(gam.*y(3)-((y(1)-ksi).^2).*y(2));
% далее запишем систем ДУ
dy(1) =(R1-2.*gam.*n.*y(1).*y(3))./R2;
dy(2) =(y(2).*(2.*y(1).*y(2).*n.*((y(1)-ksi).^2)-R1))./(R2.*(y(1)-ksi));
dy(3)=((y(2).*y(3))./R2).*(2.*(n.*gam.*y(1)+ksi.*(n-1)).*(y(1)-ksi)-(n-1).*gam.*ksi.*y(1));

six.m

Matlab M

clear all; clc;
ksi=[1 10];
y0(1)=0.75;
y0(2)=4;
y0(3)=0.75;
[T,Y]=ode23(@five,ksi,y0);
figure
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'--')
grid on
title('График искомых функций');
legend('M','Del','P');
xlabel('ksi');
disp('Вывод значений функций')
[Y(:,1),Y(:,2),Y(:,3)];

rungekutta1.m

Matlab M

%Метод Рунге-Кутта
x(1)=0;
y4(1,1)=0.75;
y4(2,1)=4;
y4(3,1)=0.75; 
b=0.5;
n=(b-x(1))/h+1;
for i=2:n
    r1=h*@five(x(i-1),y4(:,i-1));
    r2=h*@five((x(i-1)+(h/2)),(y4(:,i-1)+(r1/2)));
    r3=h*@five((x(i-1)+(h/2)),(y4(:,i-1)+(r2/2)));
    r4=h*@five((x(i-1)+h),(y4(:,i-1)+r3));
    y4(:,i)=y4(:,i-1)+((r1+(2*r2)+(2*r3)+r4)./6);
    x(i)=x(i-1)+h;
end
subplot(2,2,1);
plot(x,y4);
grid on
title('Метод Рунге-Кутта');

Только при запуске ругается на 9 строчку:

Matlab M

1 2	Error: File: rungekutta1.m Line: 9 Column: 15 Unbalanced or unexpected parenthesis or bracket.

@SSC · 15.06.2017, 11:54

Вот код в одном файле

Matlab M

function reshenie
clear all; clc;
%Метод Рунге-Кутта
x(1)=0;
y4(1,1)=0.75;
y4(2,1)=4;
y4(3,1)=0.75; 
b=0.5;
h=0.001;
n=(b-x(1))/h+1;
for i=2:n
    r1=h*five(x(i-1),y4(:,i-1));
    r2=h*five((x(i-1)+(h/2)),(y4(:,i-1)+(r1/2)));
    r3=h*five((x(i-1)+(h/2)),(y4(:,i-1)+(r2/2)));
    r4=h*five((x(i-1)+h),(y4(:,i-1)+r3));
    y4(:,i)=y4(:,i-1)+((r1+(2*r2)+(2*r3)+r4)./6);
    x(i)=x(i-1)+h;
end
plot(x,y4(1,:),x,y4(2,:),x,y4(3,:));
grid on
title('Метод Рунге-Кутта');
 
end
 
function dy=five(ksi,y)
dy = zeros(3,1);    % вектор колонка с нулевыми элементами
n=0.688377;
gam=5./3;
R1=(n-1).*ksi.*((y(1)-ksi).*y(2).*y(1)-2.*y(3));
R2=n.*ksi.*(gam.*y(3)-((y(1)-ksi).^2).*y(2));
% далее запишем систем ДУ
dy(1) =(R1-2.*gam.*n.*y(1).*y(3))./R2;
dy(2) =(y(2).*(2.*y(1).*y(2).*n.*((y(1)-ksi).^2)-R1))./(R2.*(y(1)-ksi));
dy(3)=((y(2).*y(3))./R2).*(2.*(n.*gam.*y(1)+ksi.*(n-1)).*(y(1)-ksi)-(n-1).*gam.*ksi.*y(1));
end

Однако когда стартуете при х=0, то в функции R1=0 и R2=0 и происходит деление на 0.

@Jack Frost · 15.06.2017, 12:03 **[ТС]**

Сообщение от SSC

Однако когда стартуете при х=0, то в функции R1=0 и R2=0 и происходит деление на 0.

Т.е. нужно подбирать x(1)-?

@SSC · 15.06.2017, 12:21

Сообщение от Jack Frost

Т.е. нужно подбирать x(1)-?

Почему "подбирать"?
Просто исходная система не решается в точке х=0 (посмотрите у Вас кси в знаменателе)

@Jack Frost · 15.06.2017, 12:30 **[ТС]**

Значит значения для x должны соответствовать ksi, т.е. от 1 до 10?

@SSC · 15.06.2017, 12:38

Сообщение от Jack Frost

Значит значения для x должны соответствовать ksi, т.е. от 1 до 10?

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

@Jack Frost · 22.06.2017, 20:29 **[ТС]**

Изменил значения

Matlab M

1 2	x(1)=0; b=0.5;

На

Matlab M

1 2	x(1)=1; b=10;

Вот окончательный вариант, который заработал

Matlab M

function reshenie
clear all; clc;
%Метод Рунге-Кутта
x(1)=1;
y4(1,1)=0.75;
y4(2,1)=4;
y4(3,1)=0.75; 
b=10;
h=0.001;
n=(b-x(1))/h+1;
for i=2:n
    r1=h*five(x(i-1),y4(:,i-1));
    r2=h*five((x(i-1)+(h/2)),(y4(:,i-1)+(r1/2)));
    r3=h*five((x(i-1)+(h/2)),(y4(:,i-1)+(r2/2)));
    r4=h*five((x(i-1)+h),(y4(:,i-1)+r3));
    y4(:,i)=y4(:,i-1)+((r1+(2*r2)+(2*r3)+r4)./6);
    x(i)=x(i-1)+h;
end
plot(x,y4(1,:),x,y4(2,:),x,y4(3,:));
grid on
title('Метод Рунге-Кутта');
 
end
 
function dy=five(ksi,y)
dy = zeros(3,1);    % вектор колонка с нулевыми элементами
n=0.688377;
gam=5./3;
R1=(n-1).*ksi.*((y(1)-ksi).*y(2).*y(1)-2.*y(3));
R2=n.*ksi.*(gam.*y(3)-((y(1)-ksi).^2).*y(2));
% далее запишем систем ДУ
dy(1) =(R1-2.*gam.*n.*y(1).*y(3))./R2;
dy(2) =(y(2).*(2.*y(1).*y(2).*n.*((y(1)-ksi).^2)-R1))./(R2.*(y(1)-ksi));
dy(3)=((y(2).*y(3))./R2).*(2.*(n.*gam.*y(1)+ksi.*(n-1)).*(y(1)-ksi)-(n-1).*gam.*ksi.*y(1));
end

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

M, delta, P можно взять в виде

Matlab M

1
2
3

M=y4(1,:);
delta=y4(2,:);
P=y4(3,:);

?

И как можно выделить ksi?

@SSC · 19.07.2017, 07:27

Сообщение от Jack Frost

И как можно выделить ksi?

У Вас ksi в функции five соответствуют х в головной функции

@Jack Frost · 19.07.2017, 11:24 **[ТС]**

Т.е. можно попробовать задать перед новым циклом в виде ksi=x(1); ?
А если нужно задать ksi от 1 до 10?

@SSC · 19.07.2017, 14:15

Сообщение от Jack Frost

можно попробовать задать перед новым циклом в виде ksi=x(1); ?

Я Вас не понимаю.
ksi является входной переменной в функции five определения значений производных на каждом шаге интегрирования.
Функция five вызывается при реализации численного интегрирования методом Рунге-Кутта 4-го порядка.
В какое место Вы хотите запихнуть ksi=x(1);

Сообщение от Jack Frost

А если нужно задать ksi от 1 до 10?

Ну так и задайте пределы интегрирования от 1 до 10, только при этом и начальные значения интегрируемых переменных должны быть при ksi=1, что в общем то и сделано в коде сообщения 13.

Сообщение от Jack Frost

x(1)=1;

Сообщение от Jack Frost

b=10;

@Jack Frost · 20.07.2017, 23:14 **[ТС]**

В задаче нужно реализовать циклы, один из которых:

Matlab M

tf=0.5163;
t=0.0;
t0=0.0;
dt=1.0e-2;
dksi2=1.0e-3;
i=1.0;
 
while t<0.51
    ksi2=1.0;
    yp=((tf-t0)/(tf-t))^(3*n)-1;
    while (yp)>0
        ksi2=ksi2+dksi2; 
        q=quadv(@(x)integral(ksi,delta,x),1.0,ksi2,1e-8);
        yp=((tf-t0)/(tf-t))^(3*n)-1-3*q;
    end;
    
    tout(i)=t;
    ksi2out(i)=ksi2;
    t=t+dt;
    i=i+1;
    figure(2);
    
    plot(tout,ksi2out,'r');
    xlabel('t');
    ylabel('ksi');
    pause(0.000001);
    grid on;
    end;

Matlab M

1 2	function Y = integral(ksi,delta,x) Y = interp1(ksi,delta,x,'pchip')xx;

И еще несколько циклов, где присутствуют M, P, delta, ksi.
Поэтому и задался вопросом как это лучше сделать

@Jack Frost · 24.07.2017, 14:06 **[ТС]**

Все еще актуально.

@SSC $Эксперт по математике/физике$ 3389 / 1912 / 571 Регистрация: 09.04.2015 Сообщений: 5,365
	15.06.2017, 11:19	6
	Сообщение от Jack Frost из первого файла мне не нужно, т.к. все есть в функции. Да. Сообщение от Jack Frost Задается через y4(1,1); y4(2,1); y4(3,1); ? Да. Сообщение от Jack Frost А почему именно? Определение производной происходит в предыдущей точке расчета. Обратите внимание на цикл от 2 до n 1

@Jack Frost 0 / 0 / 0 Регистрация: 26.12.2012 Сообщений: 28
	15.06.2017, 12:03 [ТС]	9
	Сообщение от SSC Однако когда стартуете при х=0, то в функции R1=0 и R2=0 и происходит деление на 0. Т.е. нужно подбирать x(1)-? 0

@SSC $Эксперт по математике/физике$ 3389 / 1912 / 571 Регистрация: 09.04.2015 Сообщений: 5,365
	15.06.2017, 12:21	10
	Сообщение от Jack Frost Т.е. нужно подбирать x(1)-? Почему "подбирать"? Просто исходная система не решается в точке х=0 (посмотрите у Вас кси в знаменателе) 1

@Jack Frost 0 / 0 / 0 Регистрация: 26.12.2012 Сообщений: 28
	15.06.2017, 12:30 [ТС]	11
	Значит значения для x должны соответствовать ksi, т.е. от 1 до 10? 0

@SSC $Эксперт по математике/физике$ 3389 / 1912 / 571 Регистрация: 09.04.2015 Сообщений: 5,365
	15.06.2017, 12:38	12
	Сообщение от Jack Frost Значит значения для x должны соответствовать ksi, т.е. от 1 до 10? Наверное да, я же не знаю физики процессов что Вы моделируете. 1

@SSC $Эксперт по математике/физике$ 3389 / 1912 / 571 Регистрация: 09.04.2015 Сообщений: 5,365
	19.07.2017, 07:27	14
	Сообщение от Jack Frost И как можно выделить ksi? У Вас ksi в функции five соответствуют х в головной функции 1

@Jack Frost 0 / 0 / 0 Регистрация: 26.12.2012 Сообщений: 28
	19.07.2017, 11:24 [ТС]	15
	Т.е. можно попробовать задать перед новым циклом в виде ksi=x(1); ? А если нужно задать ksi от 1 до 10? 0

@SSC $Эксперт по математике/физике$ 3389 / 1912 / 571 Регистрация: 09.04.2015 Сообщений: 5,365
	19.07.2017, 14:15	16
	Сообщение от Jack Frost можно попробовать задать перед новым циклом в виде ksi=x(1); ? Я Вас не понимаю. ksi является входной переменной в функции five определения значений производных на каждом шаге интегрирования. Функция five вызывается при реализации численного интегрирования методом Рунге-Кутта 4-го порядка. В какое место Вы хотите запихнуть ksi=x(1); Сообщение от Jack Frost А если нужно задать ksi от 1 до 10? Ну так и задайте пределы интегрирования от 1 до 10, только при этом и начальные значения интегрируемых переменных должны быть при ksi=1, что в общем то и сделано в коде сообщения 13. Сообщение от Jack Frost x(1)=1; Сообщение от Jack Frost b=10; 0

@Jack Frost 0 / 0 / 0 Регистрация: 26.12.2012 Сообщений: 28
	24.07.2017, 14:06 [ТС]	18
	Все еще актуально. 0