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

@S_el · 28.10.2014, 18:28

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

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

@Maikl2020 · 28.10.2014, 19:50 **[ТС]**

S_el, Да.
Думаю убрать все итерации, только метод прогонка ставит может работает.

@S_el · 28.10.2014, 19:59

Если оставить одну итерацию получится такой график:

@Maikl2020 · 28.10.2014, 20:06 **[ТС]**

S_el, А если оба итерацию убрать как получается тоже посмотрите пожалуйста если у Вас есть время.

Код добавите пожалуйста, этого графика.

@S_el · 28.10.2014, 20:11

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

А если оба итерацию убрать как получается тоже посмотрите пожалуйста если у Вас есть время.

Поясните вопрос,если не трудно.

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

Код добавите пожалуйста, этого графика.

Код я брал старый из этой темы,поэтому там может быть несоответствие некоторых коэффициентов,но основное останется:

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
86
87
88
89
90
91
92
93
function [T,x] = funcpnew(t_end)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
count=0;%[B]счетчик числа итераций[/B]
R1=0.3; 
v=1;
cp=1000;
R2=1500;
c2=660; 
L=0.2;
mf=50;
eps=1e-5;
alw=50; 
Rw=0.05; % где используется?
ro=Rw;
Tv=315;% где используется?
la2=4; 
la1=0.084; 
al1=0.5; 
al2=0.5;
et0=0.077; 
R=8.314;
E=12; 
Q=151;
k0=50;
 
T=zeros(1,mf);
Ts=zeros(1,mf);
Tn=zeros(1,mf);
alfa=zeros(1,mf);
beta=zeros(1,mf);
 
N=50;% количество пространственных узлов
T0=1200;% - начальная температура
h=L/(N-1);% шаг по x
tau= t_end/100;%шаг по t
 
Beta=(v*R1*cp)/h;
Ceta=(al1*la1+al2*la2)/(h*h);
Aeta=R2*c2+R1*cp;
Deta=-alw;%alpha0?
  for i=1: N 
    T(i) =T0;
  end;
  time =0;
while (time<t_end)
if(count>=1) % если одна итерация прошла -> прерываем выполнение цикла.
   break;
end
    count=count+1;
 time=time+tau;
alfa(1) =0;
beta(1)=300;
   for i =1 :N 
   Tn(i)=T(i);
   end;
  while(1==1)
for i=1:N  
Ts(i)=T(i);
end
  for i=2:N-1 
 
  ai =Beta-Ceta;
  bi =(Aeta/tau-Beta/h+2*Ceta/(h*h)+Deta);
  ci =Ceta/(h*h);
  fi =alw*T0-Aeta*Tn(i)/tau+ (T0*et0*Q*k0*ro*exp(-E/(R*T(i))))/T(i);
  
alfa(i) =ai/(bi-ci*alfa(i-1));
beta(i)=(ci*beta(i-1)-fi)/(bi-ci*alfa(i-1));
  end;
   T(N)=300;
         
  for i=N-1:-1:1  
T(i) =alfa(i)*T(i+1)+beta(i);
  end;
max=abs(T(1)-Ts(1));
  for i=2:N  
     if(max<abs(T(i)-Ts(i))) 
         max =abs(T(i)-Ts(i));
     end
  end
     if(max<=eps) 
         break;
     end
  
  end
end;
for i=1:N
    x(i)=h*(i-1);
end
T=T-273;
 
end

@Maikl2020 · 28.10.2014, 20:29 **[ТС]**

S_el, Думаю вот это цикл вообще не нужны

Matlab M
1
2
3
4
5
6
7
8
9
max=abs(T(1)-Ts(1));
  for i=2:N  
     if(max<abs(T(i)-Ts(i))) 
         max =abs(T(i)-Ts(i));
     end
  end
     if(max<=eps) 
         break;
     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
86
87
88
89
90
91
92
93
94
function [T,x] = funcpnew(t_end)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
count=0;%[B]счетчик числа итераций[/B]
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;
E=126000; 
Q=1513000;
k0=50000000000;
u=0.001;
v=0.5;
u0=v/u;
Al0=(2*alw)/Rw;
Fi=(R2*c2)/(R1*cp);
k=(al1*la1+al2*la2)/(R2*c2+R1*cp);
A0=Al0/(R2*c2+R1*cp);
A1=(R1*k0*Q*et0)/(R2*c2+R1*cp);
A2=(v)/(1+Fi);
%mf=500;
 
T=zeros(1,mf);
Ts=zeros(1,mf);
Tn=zeros(1,mf);
alfa=zeros(1,mf);
beta=zeros(1,mf);
 
N=100;% количество пространственных узлов
T0=300;% - начальная температура
h=L/(N-1);% шаг по x
tau= t_end/100;%шаг по t
 
Beta=(v*R1*cp)/h;
Ceta=(al1*la1+al2*la2)/(h*h);
Aeta=R2*c2+R1*cp;
Deta=-alw;%alpha0?
  for i=1: N 
    T(i) =T0;
  end;
  time =0;
while (time<t_end)
if(count>=1) % если одна итерация прошла -> прерываем выполнение цикла.
   break;
end
    count=count+1;
 time=time+tau;
alfa(1) =0;
beta(1)=1200;
   for i =1 :N 
   Tn(i)=T(i);
   end;
  for i=1:N  
Ts(i)=T(i);
end
  for i=2:N-1 
   ai =k/(h*h);
  bi =1/tau+(2*k)/(h*h)+A2/h;
  ci =A2/h+(k)/(h*h);
  fi =-1/tau+A0*(T(i)-T0)-A1*exp(-E/R/T(i));
alfa(i) =ai/(bi-ci*alfa(i-1));
beta(i)=(ci*beta(i-1)-fi)/(bi-ci*alfa(i-1));
  end;
   T(N)=300;
         
  for i=N-1:-1:1  
T(i) =alfa(i)*T(i+1)+beta(i);
  end;
max=abs(T(1)-Ts(1));
  for i=2:N  
     if(max<abs(T(i)-Ts(i))) 
         max =abs(T(i)-Ts(i));
     end
  end
     if(max<=eps) 
         break;
     end
  
  end
end;
for i=1:N
    x(i)=h*(i-1);
end
T=T-273;
end

Matlab M
1
2
3
[T,x]=funcgnew(600);
plot(x,T)
grid on

@S_el · 28.10.2014, 20:34

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

Думаю вот это цикл вообще не нужны

Я считаю,что лучше из рабочей программы никаких кусков не удалять.
Сейчас мы ищем рабочее решение,когда оно будет найдено -> удалим лишнее и займемся оптимизацией.

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

И еще у меня программа все время считает свои коэффициенты подставил

Тогда через пару часиков посмотрю.

@Maikl2020 · 29.10.2014, 12:26 **[ТС]**

S_el, Да но мой уравнения не надо проверит

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

max=abs(T(1)-Ts(1));
* for i=2:N *
* * *if(max<abs(T(i)-Ts(i)))
* * * * *max =abs(T(i)-Ts(i));
* * *end
* end
* * *if(max<=eps)
* * * * *break;
* * *end

потому, что у меня все точна задано начальные и граничное значения и константы.
По этому без эти лишний проверки программа должен работать.
Ну если у Вас будет время посмотрите пожалуйста, проста очень нуждаюсь.

Добавлено через 15 часов 47 минут
S_el, Я чуть чуть изменил параметры константы. Вот такой график получил. Но x начинается от 0
а на программа у нас 1200 почему так получается?

Кликните здесь для просмотра всего текста

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
function [T,x] = funcp1(t_end)
%Константы 
L=0.2;
mf=500;
eps=1e-5;
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;
v=1;
E=126000; 
Q=1513000;
k0=50000000000;
Al0=(2*alw)/Rw;
FFi=R2*c2/R1/cp;
kap=(al1*la1+al2*la2)/(R2*c2+R1*cp);
a0=Al0/(R2*c2+R1*cp);
a1=(R1*k0*Q*et0)/(R2*c2+R1*cp);
a2=v/(1+FFi);
T=zeros(1,mf);
Ts=zeros(1,mf);
Tn=zeros(1,mf);
alfa=zeros(1,mf);
beta=zeros(1,mf);
 
N=500;% количество пространственных узлов
T0=298;% - начальная температура
h=L/(N-1);% шаг по x
tau= t_end/100;%шаг по t
 
  for i=1: N 
    T(i) =T0;
  end;
  time =0;
while (time<t_end)
time=time+tau;
alfa(1) =0;
beta(1)=1200;
  while(1==1)
for i=1:N  
Ts(i)=T(i);
end
  for i=2:N-1 
  ai =kap/(h*h);
  bi =1/tau+(2*kap)/(h*h)+a2/h;
  ci =a2/h+kap/(h*h);
  fi =-T(i)/tau+a0*(T(i)-T0)-a1*exp(-E/(R*T(i)));
alfa(i) =ai/(bi-ci*alfa(i-1));
beta(i)=(ci*beta(i-1)-fi)/(bi-ci*alfa(i-1));
  end;
        
  for i=N-1:-1:1  
T(i) =alfa(i)*T(i+1)+beta(i);
  end;
max=abs(T(1)-Ts(1));
  for i=2:N  
     if(max<abs(T(i)-Ts(i))) 
         max =abs(T(i)-Ts(i));
     end
  end
     if(max<=eps) 
         break;
     end
  
  end
  
end;
for i=1:N
    x(i)=h*(i-1);
end
T=T-273;
 
end

Matlab M
1
2
[T,x]=funcp1(600);
>> plot(x,T)

@Maikl2020 · 29.10.2014, 12:34 **[ТС]**

S_el, График вот так получается.

@S_el · 29.10.2014, 12:36

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

Но x начинается от 0
а на программа у нас 1200 почему так получается?

Дело в масштабе,посмотрите шаг по вертикальной оси.

@Maikl2020 · 29.10.2014, 12:44 **[ТС]**

S_el, Вы не посмотрели да, чтобы убрать итерацию.

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

max=abs(T(1)-Ts(1));
* for i=2:N *
* * *if(max<abs(T(i)-Ts(i)))
* * * * *max =abs(T(i)-Ts(i));
* * *end
* end
* * *if(max<=eps)
* * * * *break;
* * *end

@S_el · 04.11.2014, 16:58

Вот,что отображается:

@Maikl2020 · 04.11.2014, 17:10 **[ТС]**

S_el, Именно мен такой нужен, отправите пожалуйста Вашу М-файл и командный код у меня так не получается.

@S_el · 04.11.2014, 17:13

вызов:

Matlab M
1
2
3
4
5
 [T,x]=func_d2(600);
T2=func_d2(1800);
T3=func_d2(3600);
plot(x,T,x,T2,x,T3)
grid on

функция:

Кликните здесь для просмотра всего текста

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
function [T,x]  = func_d2(t_end)
mf=500;
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;
T0  =300;
cp  =1000;
v  =1;
E  =126000;
Q  =1513000;
k0  =50000000000;
L  =0.1;
t1  =300;
N=500;
 
T=zeros(1,mf);
alfa=zeros(1,mf);
beta=zeros(1,mf);
 
 
Al0 =(2*alw)/rw;
FFi =r2*c2/r1/cp;
kap =(al1*la1+al2*la2)/(r2*c2+r1*cp);
a0 =Al0/(r2*c2+r1*cp);
a1 =(r1*k0*Q*et0)/(r2*c2+r1*cp);
a2 =v/(1+FFi);
alfa(1) =0;
beta(1) =1200;
 
h =L/(N-1);
tau =(t_end)/(100);
for i =1:mf
  T(i) =0;
  alfa(i) =0;
  beta(i) =0;
end
for i =1 :N 
T(i) =T0;
end
time =0;
 
for j =0 :3 
time =time+0.04;
while (time<t_end) 
    
  time =time+tau;
  alfa(1) =0;
  beta(1) =1200;
  for i =2 : N-1
    ai =kap/(h*h);
    bi =1/tau+(2*kap)/(h*h)+a2/h;
    ci =a2/h+kap/(h*h);
    fi =-T(i)/tau+a0*(T(i)-T0)-a1*exp(-E/(R*T(i)));
    alfa(i) =ai/(bi-ci*alfa(i-1));
    beta(i) =(ci*beta(i-1)-fi)/(bi-ci*alfa(i-1));
  end;
end;
 T(N) =300;
 for i = N-1:-1: 1 
 T(i) =alfa(i)*T(i+1)+beta(i);
 end;
 
 for i =1 : N
x(i)=h*(i-1);
 end;
end;
 
end

@Maikl2020 · 04.11.2014, 17:25 **[ТС]**

S_el, Вот конце концов мой тема, благодаря Вам друзья S_el, Зосима,
и сайт https://www.cyberforum.ru решился, огромный Вам спасибо.

Имена S_el, Вам огромный благодарю, с Вашей помощь это тема получил результат.

И я очень рад, что в мери есть люди как Вы.
Пусть всегда Вы помогли людям, и принести им счастья.

Новые блоги и статьи Все статьи Все блоги /
Саморегулирующийся социальный контракт для сервера 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С. Задача: Напишите приложение-калькулятор, которое помогает рассчитывать параметры кредита для аннуитетного и дифференцированного видов. . .

@S_el 2444 / 1842 / 406 Регистрация: 15.12.2013 Сообщений: 8,243
	28.10.2014, 18:28
	Maikl2020, посмотрел,проблема возникает уже на второй итерации. 1

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,451
	28.10.2014, 19:50 [ТС]
	S_el, Да. Думаю убрать все итерации, только метод прогонка ставит может работает. 0

@S_el 2444 / 1842 / 406 Регистрация: 15.12.2013 Сообщений: 8,243
	28.10.2014, 19:59
	Если оставить одну итерацию получится такой график: Миниатюры 1

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,451
	28.10.2014, 20:06 [ТС]
	S_el, А если оба итерацию убрать как получается тоже посмотрите пожалуйста если у Вас есть время. Код добавите пожалуйста, этого графика. 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,451
	29.10.2014, 12:34 [ТС]
	S_el, График вот так получается. Миниатюры 0

@S_el 2444 / 1842 / 406 Регистрация: 15.12.2013 Сообщений: 8,243
	04.11.2014, 16:58
	Вот,что отображается: Миниатюры 1

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,451
	04.11.2014, 17:10 [ТС]
	S_el, Именно мен такой нужен, отправите пожалуйста Вашу М-файл и командный код у меня так не получается. 0

@Maikl2020 Всегда онлайн 49 / 49 / 10 Регистрация: 13.04.2014 Сообщений: 1,451
	04.11.2014, 17:25 [ТС]
	S_el, Вот конце концов мой тема, благодаря Вам друзья S_el, Зосима, и сайт https://www.cyberforum.ru решился, огромный Вам спасибо. Имена S_el, Вам огромный благодарю, с Вашей помощь это тема получил результат. И я очень рад, что в мери есть люди как Вы. Пусть всегда Вы помогли людям, и принести им счастья. 0