Нелинейное дифференциальное уравнение - метод Рунге-Кутта

@Maikl2020 · Регистрация: 13.04.2014

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

Здравствуйте!
Пожалуйста подскажите тут где может быт ощипка? Вот это уравнения
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\large \dfrac{d^2T}{dx^2}=ua_{1}\dfrac{dT}{dx}+a_{2}(T-T_{0})-\dfrac{a_{3}}{T}\exp\left(-\dfrac{E}{RT}\right)$

В М-файле с именем Rk.m пишем:

Matlab M
1
2
3
4
function dT=RK(x,T)
dT=y;
dy=u*a1*T+a2*(T-T0)-(a3/T)*exp(-E/R/T);
end

Потом в командном окне вызываем функцию ode45:

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
>> alw=50; 
Rw=0.3; 
R1=0.3; 
R2=1500;
c2=660; 
la2=4; 
la1=0.084; 
al1=0.5; 
al2=0.5;
et0=0.077; 
R=8.314;
cp=1000; 
T0=300; 
E=126000; 
Q=15130000;
k0=50000000000;
u=0.0001;
v=1;
u0=v/u;
Al0=(2*Alw)/Rw;
a1=(R1*(u0+1)*cp+R2*c2*u0)/(al1*la1+al2*la2);
a2=Al0/(al1*la1+al2*la2);
a3=(R1*Q*et0*k0)/(al1*la1+al2*la2);
[T,y]=ode45(@RK,[1 10], [300, 1200]);
plot(T,y, '-k')

Зосима · 29.09.2014, 16:19

не совсем так:
файл RK.m:

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
function dT=RK(x,T)
alw=50; 
Rw=0.3; 
R1=0.3; 
R2=1500;
c2=660; 
la2=4; 
la1=0.084; 
al1=0.5; 
al2=0.5;
et0=0.077; 
R=8.314;
cp=1000; 
T0=300; 
E=126000; 
Q=15130000;
k0=50000000000;
u=0.0001;
v=1;
u0=v/u;
Al0=(2*alw)/Rw;
a1=(R1*(u0+1)*cp+R2*c2*u0)/(al1*la1+al2*la2);
a2=Al0/(al1*la1+al2*la2);
a3=(R1*Q*et0*k0)/(al1*la1+al2*la2);
 
dT(1,1) = T(2); % первая производная dT/dx
dT(2,1)=u*a1*T(2) + a2*(T(1)-T0) - (a3/T(1))*exp(-E/R/T(1));
end

и запуск расчетов:

Matlab M
1
2
[T,y]=ode45(@RK,[1 10], [300, 1200]);
plot(T,y)

но оно почему-то улетает в бесконечность... проверь условия.

@Maikl2020 · 29.09.2014, 18:57 **[ТС]**

Зосима, Спасибо что ответили но сам тоже не понимаю почему улетает в бесконечность

Какие условия Вы имейте виду?

И еще

[T,y]
ans =
1.0e+302 *
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
Эти три столбцы какой значение он выводит?

Зосима · 29.09.2014, 19:17

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

Какие условия Вы имейте виду?

имею ввиду правильность всех формул и чисел

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

Эти три столбцы какой значение он выводит?

не знаю, так не видно

@Maikl2020 · 29.09.2014, 19:27 **[ТС]**

Зосима, Может это три столбец dT(1,1) и dT(2,1) и x да?

Зосима · 29.09.2014, 19:34

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

Может это три столбец dT(1,1) и dT(2,1) и x да?

да

только первый столбец - это значения х, второй - значения Т и третий значения dT/dx

@Maikl2020 · 30.09.2014, 07:25 **[ТС]**

Зосима, Если так то почему значения x всегда ноль а там дано же интервал, в том числе T тоже

@S_el · 30.09.2014, 08:40

Maikl2020, попробуйте решить своё уравнение явной схемой,такие решения более устойчивы.

Зосима · 30.09.2014, 09:36

Maikl2020, видишь вначале стоит множитель "1.0e+302" это 10³⁰². х таки да меняется от 0 до 10, но этого не видно, т.к. из-за множителя значущие цифры теряются

вот например: x = 10;
10
1 * 10e+1
0.1 * 10e+2
0.001 * 10e+3
0.0001 * 10e+4
0.0000 * 10e+5
Получаем нули уже при множителе 10⁵ при формате отображения 4 знака после запятой

@Maikl2020 · 30.09.2014, 16:28 **[ТС]**

S_el, Например какое методы Вы имейте виду?
например прогонка годится?

@S_el · 30.09.2014, 16:32

Maikl2020, я поторопился с выводами,не посмотрел на коэффициенты.
Лучше сделайте как сказал Зосима: перепроверьте числа.

Зосима · 30.09.2014, 17:25

Maikl2020, S_el намекал решить дифуру вручную и получить аналитическое решение(формулу)

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

Добавлено через 1 минуту
*единицы измерения в СИ перевел?

@S_el · 30.09.2014, 17:45

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

(хотя радт эексперимента можешь попробовать другие решатели, вмсто ode45 взять ode115, ode26 или какие там еще есть, не помню)

Я пробовал -> с такими коэффициентами это не важно.

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

S_el намекал решить дифуру вручную и получить аналитическое решение(формулу)

Может оно и существует,но с такими числами его не получить

Даже если отбросить член с экспонентой,деленной на y ,рассмотрев его позже.

@Maikl2020 · 01.10.2014, 19:27 **[ТС]**

S_el, Ну что у Вас получился ?

Зосима · 01.10.2014, 20:22

Maikl2020, мы от тебя ждем результатов проверки всех чисел и формул!!!!!!

@Maikl2020 · 01.10.2014, 20:28 **[ТС]**

Зосима, Все числи и формул правильно, но думаю граничный условия не так, потому, что в Delphi написал тоже так получился

@S_el · 01.10.2014, 20:32

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

Ну что у Вас получился ?

Ничего.

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

Зосима · 01.10.2014, 20:48

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

но думаю граничный условия не так

стоп, т.е. это: ode45(@RK,[1 10], [300, 1200]) не начальные, а граничные условия? (в начале и конце)

@Maikl2020 · 02.10.2014, 11:33 **[ТС]**

Зосима, S_el, правда сам путаю но график должен быт примерно такой

Нелинейное дифференциальное уравнение - метод Рунге-Кутта

@Maikl2020 · 03.10.2014, 12:11 **[ТС]**

Здравствуй Зосима, посмотри пожалуйста я чуть чуть по менял уравнения т.е. обезразмероваль.
Помоги пожалуйста код программа правильно или нет, а то ошибка дает.
Вот это уравнения
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\large \dfrac{d^2\theta}{dx^2}=L\cdot u\cdot a_{1}\dfrac{dT}{dx}+L^2*a_{2}\theta-\dfrac{a_{3}\cdot L^2}{T0+\theta(T_e-T_0)}\exp\left(-\dfrac{E}{R(T0+\theta(T_e-T_0))}\right)$
потом из уравнения получил следующий система
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\large \dfrac{d\theta}{dx}=y,\\\dfrac{dy}{dx}=L\cdot u\cdot a_{1}\cdot y+L^2\cdot a_{2}\cdot\theta-\dfrac{a_{3}\cdot L^2}{T_0+\theta(T_e-T_0)}\exp\left(-\dfrac{E}{R(T_0+\theta(T_e-T_0))}\right)$
граничный условия для система
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\large y=0.1, \\\theta=\dfrac{T_v-T_0}{T_e-T0}$

M-файл

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 dTETA=RK(x,TETA)
dTETA=zeros(2,1);
alw=50; 
Rw=0.05; 
R1=0.3; 
R2=1500;
c2=660; 
la2=4; 
la1=0.084; 
al1=0.5; 
al2=0.5;
et0=0.077; 
R=8.314;
cp=1000; 
Tv=315;
T0=300; 
E=12; 
Q=1513000;
k0=50000000000;
u=0.001;
v=0.5;
u0=v/u;
L=0.5;
Te=T0+(Q*et0)/(1+(R2*c2*u)/(R1*cp*(v+u)));
TETA=(Tv-T0)/(Te-T0);
Al0=(2*alw)/Rw;
a1=(R1*(u0+1)*cp+R2*c2*u0)/(al1*la1+al2*la2);
a2=Al0/(al1*la1+al2*la2);
a3=(R1*Q*et0*k0)/(al1*la1+al2*la2);
dTETA=Y; % первая производная dT/dx
dY=L*u*a1*y+L*L*a2*TETA-((a3*L*L)/(T0+TETA*(Te-T0)))*exp(-E/(R*(T0+TETA*(Te-T0))));
end

код

Matlab M
1
2
[x,TETA]=ode45(@RKU,[1 10], [TETA 0.1]);
plot(x,TETA,'-k')

Новые блоги и статьи Все статьи Все блоги /
YAFU@home — распределённые вычисления для математики. На CPU Programma_Boinc 20.01.2026 YAFU@home — распределённые вычисления для математики. На CPU YAFU@home — это BOINC-проект, который занимается факторизацией больших чисел и исследованием aliquot-последовательностей. Звучит. . .	http://iceja.net/ математические сервисы iceja 20.01.2026 Обновила свой сайт http:/ / iceja. net/ , приделала Fast Fourier Transform экстраполяцию сигналов. Однако предсказывает далеко не каждый сигнал (см ограничения http:/ / iceja. net/ fourier/ docs ). Также. . .	http://iceja.net/ сервер решения полиномов iceja 18.01.2026 Выкатила http:/ / iceja. net/ сервер решения полиномов (находит действительные корни полиномов методом Штурма). На сайте документация по API, но скажу прямо VPS слабенький и 200 000 полиномов. . .	Расчёт переходных процессов в цепи постоянного тока igorrr37 16.01.2026 / * Дана цепь постоянного тока с R, L, C, k(ключ), U, E, J. Программа составляет систему уравнений по 1 и 2 законам Кирхгофа, решает её и находит: токи, напряжения и их 1 и 2 производные при t = 0;. . .
Восстановить юзерскрипты Greasemonkey из бэкапа браузера damix 15.01.2026 Если восстановить из бэкапа профиль Firefox после переустановки винды, то список юзерскриптов в Greasemonkey будет пустым. Но восстановить их можно так. Для этого понадобится консольная утилита. . .	Сукцессия микоризы: основная теория в виде двух уравнений. anaschu 11.01.2026 https:/ / rutube. ru/ video/ 7a537f578d808e67a3c6fd818a44a5c4/	WordPad для Windows 11 Jel 10.01.2026 WordPad для Windows 11 — это приложение, которое восстанавливает классический текстовый редактор WordPad в операционной системе Windows 11. После того как Microsoft исключила WordPad из. . .	Classic Notepad for Windows 11 Jel 10.01.2026 Old Classic Notepad for Windows 11 Приложение для Windows 11, позволяющее пользователям вернуть классическую версию текстового редактора «Блокнот» из Windows 10. Программа предоставляет более. . .

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,445
	29.09.2014, 18:57 [ТС]
	Зосима, Спасибо что ответили но сам тоже не понимаю почему улетает в бесконечность Какие условия Вы имейте виду? И еще [T,y] ans = 1.0e+302 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Эти три столбцы какой значение он выводит? 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,445
	29.09.2014, 19:27 [ТС]
	Зосима, Может это три столбец dT(1,1) и dT(2,1) и x да? 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,445
	30.09.2014, 07:25 [ТС]
	Зосима, Если так то почему значения x всегда ноль а там дано же интервал, в том числе T тоже 0

@S_el 2444 / 1842 / 406 Регистрация: 15.12.2013 Сообщений: 8,243
	30.09.2014, 08:40
	Maikl2020, попробуйте решить своё уравнение явной схемой,такие решения более устойчивы. 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	30.09.2014, 09:36
	Maikl2020, видишь вначале стоит множитель "1.0e+302" это 10³⁰². х таки да меняется от 0 до 10, но этого не видно, т.к. из-за множителя значущие цифры теряются вот например: x = 10; 10 1 * 10e+1 0.1 * 10e+2 0.001 * 10e+3 0.0001 * 10e+4 0.0000 * 10e+5 Получаем нули уже при множителе 10⁵ при формате отображения 4 знака после запятой 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,445
	30.09.2014, 16:28 [ТС]
	S_el, Например какое методы Вы имейте виду? например прогонка годится? 0

@S_el 2444 / 1842 / 406 Регистрация: 15.12.2013 Сообщений: 8,243
	30.09.2014, 16:32
	Maikl2020, я поторопился с выводами,не посмотрел на коэффициенты. Лучше сделайте как сказал Зосима: перепроверьте числа. 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	30.09.2014, 17:25
	Maikl2020, S_el намекал решить дифуру вручную и получить аналитическое решение(формулу) однако я не уверен что такое явное решение существует... какой бы ты метод не взял: прогонки, Эйлера, Рунге-Кутта - при этих данных результат будет тоже. (хотя радт эексперимента можешь попробовать другие решатели, вмсто ode45 взять ode115, ode26 или какие там еще есть, не помню) Добавлено через 1 минуту *единицы измерения в СИ перевел? 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,445
	01.10.2014, 19:27 [ТС]
	S_el, Ну что у Вас получился ? 0

Зосима 5245 / 3573 / 379 Регистрация: 02.04.2012 Сообщений: 6,477 Записей в блоге: 18
	01.10.2014, 20:22
	Maikl2020, мы от тебя ждем результатов проверки всех чисел и формул!!!!!! 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,445
	01.10.2014, 20:28 [ТС]
	Зосима, Все числи и формул правильно, но думаю граничный условия не так, потому, что в Delphi написал тоже так получился 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,445
	02.10.2014, 11:33 [ТС]
	Зосима, S_el, правда сам путаю но график должен быт примерно такой 0