Метод Рунге-Кутты для решения системы ОДУ

@Svetl123 · Регистрация: 30.01.2012

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

4 дня сижу и пытаюсь понять этот метод(

Зосима · 10.12.2013, 20:51

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

Но вот видео с этой анимации не пойму как делать

в этом сообщении: https://www.cyberforum.ru/post3985757.html
под спойлером есть команда, которая создает видео

@Svetl123 · 10.12.2013, 23:04 **[ТС]**

Сообщение от Зосима

в этом сообщении: https://www.cyberforum.ru/post3985757.html
под спойлером есть команда, которая создает видео

Не додумался как через movie2avi, сделал через videowriter

Мне нужно сделать таблицу с погрешностями между этим решением и ode45
Посоветуешь может какую ф-ю для облегчения создания таблицы?
Или придется все в ручную?

Зосима · 11.12.2013, 13:28

Табличку нарисовать сложнее чем сделать видео

но вся загвоздка в том, что в ode шаг выбирается автоматически и он не постоянный, поэтому количество точек ручного и автоматического решений будет разным....
А!

не, массив времени можно задать в ode:

Matlab M
1
2
3
4
5
6
7
8
clear, clc
 
y0 = [0 0.1 0.1];
[T,U] = three('one',y0,0,10,80);
[t v] = ode45('one',T,y0);
d = abs(U-v);
disp('       T        X         Y         Z         dX        dY        dZ')
disp([T U d])

@Svetl123 · 11.12.2013, 15:20 **[ТС]**

Спасибо, друГ) Ноя уже сделал так с утра додумался) хы)
у меня вопросец я переработал код в апостериорный но с вайлом мучаюсь не continue не break не берет) как быть не знаю)
Два часа назад защитил работу )

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
function [T U]=two(f,y0,a,b,M)
% f - имя функции системы
% y0 - начальные условия
% a,b - интервал времени
% M - кол-во точек разбиения
h = (b-a)/M;
U = zeros(M+1,length(y0));
T = (a:h:b)';
U(1,:) = y0;
eps = 0.01;
for j=1:M 
    while(1)
    k1=h*feval(f,T(j), U(j,:) )';
    k2=h*feval(f,T(j)+h/2, U(j,:)+k1/2)';
    k3=h*feval(f,T(j)+h/2, U(j,:)+k2/2)';
    k4=h*feval(f,T(j)+h, U(j,:)+k3)';
    D = U(j,:) + (k1+2*k2+2*k3+k4)/6 ;
    h=h/2;
    k1=h*feval(f,T(j), U(j,:) )';
    k2=h*feval(f,T(j)+h/2, U(j,:)+k1/2)';
    k3=h*feval(f,T(j)+h/2, U(j,:)+k2/2)';
    k4=h*feval(f,T(j)+h, U(j,:)+k3)';
    D1 = U(j,:) + (k1+2*k2+2*k3+k4)/6 ;
    if max(abs(D-D1)) > eps;
        h=h/2;
    continue
    else
        T(j+1) = T(j) + h;
        U(j+1,:) = D;
    end
    end
    h = 2*h;
end
end

Зосима · 11.12.2013, 16:03

Svetl123, ну а где ты в своей програмке задаешь точность eps ?

Добавлено через 5 минут
ой

не заметил...

Добавлено через 10 минут
Вот, причесал:

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
function [T U]=two(f,y0,a,b,eps)
% f - имя функции системы
% y0 - начальные условия
% a,b - интервал времени
% eps - точность
h = (b-a)/50; % начальное приближение шага
T(1) = a;
U(1,:) = y0(:)' ;
j = 1;
while T(end)<=b
    k1=h*feval(f,T(j), U(j,:) )';
    k2=h*feval(f,T(j)+h/2, U(j,:)+k1/2)';
    k3=h*feval(f,T(j)+h/2, U(j,:)+k2/2)';
    k4=h*feval(f,T(j)+h, U(j,:)+k3)';
    D = U(j,:) + (k1+2*k2+2*k3+k4)/6 ;
    h=h/2; 
    k1=h*feval(f,T(j), U(j,:) )';
    k2=h*feval(f,T(j)+h/2, U(j,:)+k1/2)';
    k3=h*feval(f,T(j)+h/2, U(j,:)+k2/2)';
    k4=h*feval(f,T(j)+h, U(j,:)+k3)';
    D1 = U(j,:) + (k1+2*k2+2*k3+k4)/6 ;
    if max(abs(D-D1)) > eps; % если точность не достигнута
        h=h/2; % уменьшаем шаг
        continue
    else
        T(j+1) = T(j) + h;
        U(j+1,:) = D;
        h = 2*h;
        j = j+1;
    end
end
end

и вызов

Matlab M
1
2
3
4
5
6
7
8
clear, clc
 
y0 = [0 0.1 0.1];
[T,U] = three('one',y0,0,10,80);
[t v] = two('one',y0,0,10,0.01);
plot3(U(:,1), U(:,2), U(:,3),'r',...
    v(:,1),v(:,2),v(:,3),'b')
grid on

@Svetl123 · 11.12.2013, 16:34 **[ТС]**

Matlab M
1
2
3
4
5
6
7
8
clear, clc
 
y0 = [0 0.1 0.1];
[T,U] = three('one',y0,0,10,80);
[t v] = two('one',y0,0,10,0.01);
plot3(U(:,1), U(:,2), U(:,3),'r',...
    v(:,1),v(:,2),v(:,3),'b')
grid on

файла-ф-ии three(это у меня запускающий файл так называется главный т.е) у мну точно нет а если делать так то графики не содятся:

Matlab M
1
2
3
4
5
clear, clc
[T,U] = two('one',[0 0.1 0.1],0,10,0.01);
[t v] = ode45('one',[0 10],[0 0.1 0.1]);
plot3(U(:,1), U(:,2), U(:,3),'r',v(:,1), v(:,2), v(:,3),'bo')
grid on

@Svetl123 · 11.12.2013, 17:39 **[ТС]**

Хотя все я разобрался)
Спасибо, тебе)
Меня Никита кстати звать, на будущее, надеюсь увидимся)
В вк все ранво не отвечаешь)

Зосима · 12.12.2013, 18:29

Svetl123, молоток!

Не по теме:

в вк я заглядываю часто, но не всегда отвечаю :pardon:

Новые блоги и статьи Все статьи Все блоги /
Саморегулирующийся социальный контракт для сервера cross-section. Hrethgir 14.08.2026 С кодом конечно таких глубоких размышлений пока не было, впрочем я уже привык к алгоритмизации. Суть предмета записи: снова в диалоге с нейросетью (я взял пока себе ник для учётки админа - Rector). . . .	Часы электронные Uhbif79 12.08.2026 Выкладываю программу часов. Программа позволяет: 1. Использовать системное время и дату, 2. Есть возможность вводить время и дату вручную. 3. Реализованы 2 будильника: начало и конец рабочего дня. . . .	Часы с будильником на основе класса QLCDNumber Uhbif79 12.08.2026 Всем добрый день, выкладываю программу часов с будильником на основе класса QLCDNumber. Здесь я пробовал самостоятельно создавал классы, впервые столкнулся с видимостью переменной одного класса из. . .	Установка MinGW GCC 16.2 и CMake 8Observer8 10.08.2026 VK Видео: https:/ / vkvideo. ru/ video-240781534_456239017 YouTube: eY5-5PyI9NM Текстовая версия
Неделя из жизни имитационной модели склада: мои кривые руки растут, откуда надо anaschu 10.08.2026 Неделя из жизни имитационной модели склада: как я почти написал неправильную логику и что с этим делать Работаю сейчас над учебно-рабочим проектом: строю в AnyLogic имитационную модель процессов. . .	Калькулятор для расчета родства russiannick 07.08.2026 1. Задача: Создать калькулятор для расчета родства. Родственных связей существует 8 ступеней, такие как: p - отец P - мать q - муж Q - жена b - брат B - сестра s - сын S - дочь	Мир по моей воле kumehtar 07.08.2026 Когда-то кажется, что всё просто. Ты весь такой светлый. Причиняешь добро. Борешься за справедливость в этом тёмном мире. Потом начинаешь замечать одну неприятную вещь. Почти каждый хороший. . .	Кредитный калькулятор Maks 05.08.2026 Решение задачи по прикладной информатике средствами 1С. Задача: Напишите приложение-калькулятор, которое помогает рассчитывать параметры кредита для аннуитетного и дифференцированного видов. . .

@Svetl123 0 / 0 / 0 Регистрация: 30.01.2012 Сообщений: 16
	11.12.2013, 17:39 [ТС]
	Хотя все я разобрался) Спасибо, тебе) Меня Никита кстати звать, на будущее, надеюсь увидимся) В вк все ранво не отвечаешь) 0

Зосима 5246 / 3574 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	12.12.2013, 18:29
	Svetl123, молоток! Не по теме: в вк я заглядываю часто, но не всегда отвечаю :pardon: 0