Форум программистов, компьютерный форум, киберфорум
Matlab
Войти
Регистрация
Восстановить пароль
 
Рейтинг 4.86/7: Рейтинг темы: голосов - 7, средняя оценка - 4.86
0 / 0 / 0
Регистрация: 11.02.2014
Сообщений: 9
1

Решение системы ОДУ

11.02.2014, 22:44. Показов 1442. Ответов 5
Метки нет (Все метки)

Задача: необходимо решить систему ОДУ:
Решение системы ОДУ
,
где
Решение системы ОДУ

начальные условия
Решение системы ОДУ


Все переменные тензора Т и начальных условий расписаны в прилагаемых файлах. Сложность в том, что элементы тензора зависят от переменной (от высоты = r) по которой происходит дифференцирование, поэтому на каждом шаге вычисления тензор Т должен меняться в зависимости от высоты. Элементы тензора зависят от высоты потому что в них содержаться частота соударений и концентрация которые зависят от высоты. Для частоты и концентрации есть экспериментальные данные с шагом 2000 метров, которые описывают их зависимость от высоты, но эти данные нужно проинтерполировать, чтобы решать систему ОДУ с более высокой точностью. Диапазон изменения r от 100000 до 50000 метров, при этом решать нужно именно сверху вниз (от большего r к меньшему). Начальные условия задаются на верхней границе.

Прилагаю файлы своей программы. Программа работает, но результаты получаются неверными. Программа.rar
Вопросы: 1. может не верно реализовано изменение тензора Т в зависимости от высоты и это делается как-то иначе?
2. может задачи такого рода решаются вообще как-то иначе и даже не в матлабе?

Заранее благодарен!
0
Вложения
Тип файла: rar Вложения(система ОДУ и тд).rar (46.5 Кб, 7 просмотров)
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
11.02.2014, 22:44
Ответы с готовыми решениями:

Решение системы ОДУ
Добрый день. Извините за мой, достаточно наивный, вопрос, но как решить данную задачу, используя...

Решение системы ОДУ
Добрый день. Разбираюсь с . Нужно решить систему уравнений. Все возможные методические...

Решение системы ОДУ
Ребят, не доходит до меня способ решения, надеюсь на помощь!

Решение системы ОДУ первого порядка
Необходимо решить систему(3.11): Для этого создаю m-файл Angle.m: function dy=Angle(t,y)...

__________________
5
5142 / 3480 / 356
Регистрация: 02.04.2012
Сообщений: 6,387
Записей в блоге: 16
12.02.2014, 12:10 2
Хе, задачка не из простых, но я чуток причесал
1. функция (а не извращение ) расчета концентрации:
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
function Ni=koncentr(ri)
%Программа интерполяции концентрации
r = 50000:2000:540000; %диапазон высот для данных которые нужно проинтерполировать
N =[8000000; 10000000;12000000;16000000;20000000;31000000;43000000;60000000;80000000;98000000;117000000;200000000;400000000;600000000;1000000000;
1800000000;3000000000;5000000000;12000000000;27000000000;65000000000;80000000000;90000000000;105000000000;110000000000;120000000000;125000000000;
127000000000;130000000000;135000000000;140000000000;145000000000;147000000000;150000000000;153000000000;155000000000;154000000000;153900000000;
153000000000;152000000000;150000000000;149000000000;148000000000;147000000000;146000000000;145000000000;146000000000;147000000000;148000000000;
149000000000;150000000000;152000000000;154000000000;157000000000;158000000000;160000000000;162000000000;165000000000;167000000000;171000000000;
175000000000;177000000000;180000000000;183000000000;187000000000;190000000000;197000000000;200000000000;205000000000;210000000000;218000000000;
220000000000;225000000000;230000000000;235000000000;240000000000;242000000000;250000000000;252000000000;257000000000;260000000000;270000000000;
275000000000;280000000000;285000000000;290000000000;291000000000;300000000000;302000000000;307000000000;310000000000;319000000000;322000000000;
330000000000;332000000000;340000000000;341000000000;345000000000;347000000000;349000000000;350000000000;349700000000;349600000000;349500000000;
349400000000;349300000000;349000000000;348000000000;345000000000;343000000000;340000000000;339000000000;338000000000;335000000000;333000000000;
330000000000;327000000000;325000000000;320000000000;317000000000;315000000000;310000000000;307000000000;305000000000;300000000000;295000000000;
293000000000;290000000000;285000000000;280000000000;275000000000;270000000000;267000000000;262000000000;260000000000;255000000000;253000000000;
250000000000;245000000000;240000000000;235000000000;233000000000;230000000000;225000000000;220000000000;219000000000;215000000000;210000000000;
205000000000;203000000000;200000000000;197000000000;195000000000;190000000000;187000000000;185000000000;183000000000;180000000000;178000000000;
176000000000;175000000000;170000000000;167000000000;165000000000;163000000000;160000000000;159000000000;157000000000;155000000000;153000000000;
150000000000;147000000000;145000000000;144000000000;143000000000;140000000000;137000000000;135000000000;134000000000;132000000000;130000000000;
128000000000;126000000000;125000000000;124000000000;123000000000;122000000000;120000000000;117000000000;116000000000;115000000000;114000000000;
113000000000;112000000000;110000000000;107000000000;106000000000;105000000000;104000000000;103000000000;102000000000;100000000000;99000000000;
98000000000;96000000000;95000000000;94000000000;93000000000;91000000000;90000000000;89000000000;88000000000;87000000000;86000000000;85000000000;
84000000000;83000000000;81000000000;80000000000;79000000000;78000000000;77000000000;76000000000;75000000000;74000000000;73000000000;72500000000;
72000000000;71000000000;70000000000;69000000000;68000000000;67500000000;67000000000;66300000000;66000000000;65000000000;64000000000;63700000000;
63000000000;62000000000;61800000000;61000000000;60500000000;60000000000;59000000000];
Ni = spline(r, N, ri); % сплайны красивше
2. функция интерополяции частоты
Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function vi = chastota(ri)
%Программа интерполяции частоты
r = 50000:2000:540000; %диапазон высот для данных которые нужно проинтерполировать
v =[117000000.0;90000000.0;70000000.0;50000000.0;40000000.0;33500000.0;23000000.0;18000000.0;14000000.0;11000000.0;8530000.0;6000000.0;4000000.0;
3000000.0;2000000.0;1700000.0;980000.0;800000.0;600000.0;400000.0;297000.0;230000.0;180000.0;130000.0;100000.0;79000.0;56000.0;48000.0;35000.0;
25000.0;20000.0;17500.0;12000.0;8000.0;6400.0;5320.0;4800.0;3800.0;3300.0;2900.0;2500.0;2200.0;2000.0;1800.0;1600.0;1450.0;1400.0;1300.0;1200.0;
1130.0;1070.0;1030.0;960.0;920.0;890.0;877.0;840.0;820.0;790.0;770.0;750.0;730.0;710.0;690.0;680.0;667.0;658.0;659.0;660.0;667.0;670.0;677.0;
680.0;690.0;695.0;696.0;690.0;685.0;680.0;670.0;660.0;650.0;645.0;640.0;620.0;600.0;590.0;580.0;570.0;560.0;550.0;530.0;520.0;510.0;500.0;490.0;
475.0;465.0;455.0;445.0;435.0;420.0;410.0;390.0;380.0;370.0;360.0;350.0;340.0;330.0;320.0;310.0;300.0;290.0;280.0;270.0;260.0;255.0;250.0;240.0;
230.0;225.0;220.0;210.0;205.0;196.0;193.0;190.0;185.0;180.0;175.0;173.0;170.0;165.0;160.0;155.0;153.0;150.0;147.0;143.0;140.0;137.0;135.0;130.0;
127.0;125.0;123.0;120.0;117.0;115.0;114.0;107.0;105.0;103.0;101.0;100.0;98.0;96.0;94.0;92.0;91.0;89.0;88.0;86.0;85.0;83.0;81.0;80.0;79.0;78.0;
76.0;75.0;74.0;73.0;72.0;70.0;69.0;68.0;67.5;66.5;66.0;65.0;64.5;64.0;63.5;63.0;62.0;61.0;60.0;59.0;58.0;57.0;56.5;56.0;55.7;55.0;54.0;53.0;52.5;
52.0;51.0;50.0;49.8;49.0;48.0;47.0;46.5;46.0;45.5;45.0;44.0;43.5;43.0;42.0;41.5;41.0;40.5;40.0;39.0;38.5;38.0;37.5;37.0;36.2;36.0;35.0;34.7;34.5;
34.0;33.5;33.0;32.5;32.0;31.8;31.5;31.0;30.5;30.0;29.7;29.5;29.0;28.5;28.3;28.0;27.5;27.0];
vi = spline(r, v, ri);
3. функции систем:
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
function dx=uravneniyax(rd, x)
global  vni w k Y m n l
 
Ni=koncentr(rd);
vi=chastota(rd);
 
 
%Элементы матрицы
wpl=56.41460232*sqrt(Ni);
Z=vi/w;
X=(wpl^2)/(w^2);
U=1+1i*Z;
S=vni/(k*(rd+6370000));
 
%Опишем компоненты тензора диэлектрической проницаемости
E=U*((U^2)-(Y^2));
Err=((U^2)*(U-X)-(U-(n^2)*X)*(Y^2))/E;
Ert=((l*n*Y-1i*m*U)*X*Y)/E;
Erf=((m*n*Y+1i*l*U)*X*Y)/E;
 
Efr=((m*n*Y-1i*l*U)*X*Y)/E;
Eft=((l*m*Y+1i*n*U)*X*Y)/E;
Eff=((U^2)*(U-X)-(U-(m^2)*X)*(Y^2))/E;
 
Etr=((l*n*Y+1i*m*U)*X*Y)/E;
Ett=((U^2)*(U-X)-(U-(l^2)*X)*(Y^2))/E;
Etf=((l*m*Y-1i*n*U)*X*Y)/E;
 
%Опишем переменные сигма
btt=Err*Ett-Etr*Ert;
btf=Err*Etf-Etr*Erf;
bft=Err*Eft-Efr*Ert;
bff=Err*Eff-Efr*Erf;
 
%Матрица Т
T(1,1) = 1i*k*(-1)*S*Ert/Err;
T(1,2) = 1i*k*(-1)*S*Erf/Err;
T(1,3) = 0;
T(1,4) = 1i*k*(1-S^2/Err);
 
T(2,1) = 0;
T(2,2) = 0;
T(2,3) = 1i*k*(-1);
T(2,4) = 0;
 
T(3,1) = 1i*k*(-1)*bft/Err;
T(3,2) = 1i*k*(-bff/Err+S^2);
T(3,3) = 0;
T(3,4) = 1i*k*S*Efr/Err;
 
T(4,1) = 1i*k*btt/Err;
T(4,2) = 1i*k*btf/Err;
T(4,3) = 0;
T(4,4) = 1i*k*(-1)*S*Etr/Err;
 
dx = T*x;
% %Функция
% dx(1)=T11*x(1)+T12*x(2)+T13*x(3)+T14*x(4); 
% dx(2)=T21*x(1)+T22*x(2)+T23*x(3)+T24*x(4); 
% dx(3)=T31*x(1)+T32*x(2)+T33*x(3)+T34*x(4); 
% dx(4)=T41*x(1)+T42*x(2)+T43*x(3)+T44*x(4);
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
function dy=uravneniyay(rd, y)
global  vni w k Y m n l
 
Ni=koncentr(rd);
vi=chastota(rd);
 
%Элементы матрицы
wpl=56.41460232*sqrt(Ni);
Z=vi/w;
X=(wpl^2)/(w^2);
U=1+1i*Z;
S=vni/(k*(rd+6370000));
 
%Опишем компоненты тензора диэлектрической проницаемости
E=U*((U^2)-(Y^2));
Err=((U^2)*(U-X)-(U-(n^2)*X)*(Y^2))/E;
Ert=((l*n*Y-1i*m*U)*X*Y)/E;
Erf=((m*n*Y+1i*l*U)*X*Y)/E;
 
Efr=((m*n*Y-1i*l*U)*X*Y)/E;
Eft=((l*m*Y+1i*n*U)*X*Y)/E;
Eff=((U^2)*(U-X)-(U-(m^2)*X)*(Y^2))/E;
 
Etr=((l*n*Y+1i*m*U)*X*Y)/E;
Ett=((U^2)*(U-X)-(U-(l^2)*X)*(Y^2))/E;
Etf=((l*m*Y-1i*n*U)*X*Y)/E;
 
%Опишем переменные сигма
btt=Err*Ett-Etr*Ert;
btf=Err*Etf-Etr*Erf;
bft=Err*Eft-Efr*Ert;
bff=Err*Eff-Efr*Erf;
 
%Матрица Т
T(1,1) = 1i*k*(-1)*S*Ert/Err;
T(1,2) = 1i*k*(-1)*S*Erf/Err;
T(1,3) = 0;
T(1,4) = 1i*k*(1-S^2/Err);
 
T(2,1) = 0;
T(2,2) = 0;
T(2,3) = 1i*k*(-1);
T(2,4) = 0;
 
T(3,1) = 1i*k*(-1)*bft/Err;
T(3,2) = 1i*k*(-bff/Err+S^2);
T(3,3) = 0;
T(3,4) = 1i*k*S*Efr/Err;
 
T(4,1) = 1i*k*btt/Err;
T(4,2) = 1i*k*btf/Err;
T(4,3) = 0;
T(4,4) = 1i*k*(-1)*S*Etr/Err;
 
dy = T*y;
% %Функция
% dy(1)=T11*y(1)+T12*y(2)+T13*y(3)+T14*y(4); 
% dy(2)=T21*y(1)+T22*y(2)+T23*y(3)+T24*y(4); 
% dy(3)=T31*y(1)+T32*y(2)+T33*y(3)+T34*y(4); 
% dy(4)=T41*y(1)+T42*y(2)+T43*y(3)+T44*y(4);
 
end
4. константы ввел в тело программы расчета:
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
clear; clc; % Очистка памяти и экрана
format long
global  vni w k Y m n l
 
f=500; %рабочая частота
w=2*pi*f; %частота
Y=11097286.68/w; 
Yr=-10661904.18/w; %для продольной Y 
vni=77.42028+1i*6.982384;
k=w/299792458;
 
l=0.277357634;
m=0;
n=-0.960766847; 
%Концентрация и частота соударений для начальных условий
NN=1.2*10^11;
VV=79000;
XX=3182.607355*NN/(w^2);
n0=(1-(XX/(1+(1i*VV/w)+Yr)))^0.5;
nx=(1-(XX/(1+(1i*VV/w)-Yr)))^0.5;
 
h =100; % Шаг интегрирования
min=50000; %начальная высота
max=100000;%конечная высота
 
 
ri =max:-h:min; %диапазон интегрирования
 
%Начальные условия------------0-------------------------------------
x0=[1/10^0; -1i/10^0; 1i*n0/10^0; n0/10^0]; %Начальные значения
 
xn=[1/10^0; 1i/10^0; -1i*nx/10^0; nx/10^0]; %Начальные значения
% ------------------------------------------------------------------
 
options=odeset('AbsTol',1e-30);
 
[rd,X]=ode45('uravneniyax',ri,x0,options);
[rd,Y]=ode45('uravneniyay',ri,xn,options);
 
 
% Графики-----------------------------------------------------------
% %1 
%  semilogy(rd,abs(X(:,1)),'b',rd,abs(Y(:,1)),'r','LineWidth',2); 
%  title('et') 
%  xlabel('r') 
%  ylabel('|et|') 
%  legend('eto','ete')
% 
% %2
% figure;
%  semilogy(rd,abs(X(:,2)),'b',rd,abs(Y(:,2)),'r','LineWidth',2); 
%  title('ef') 
%  xlabel('r') 
%  ylabel('|ef|') 
%  legend('efo','efe')
%  
%  %3
% figure;
%  semilogy(rd,abs(X(:,3)),'b',rd,abs(Y(:,3)),'r','LineWidth',2);
%  title('ht') 
%  xlabel('r') 
%  ylabel('|ht|') 
%  legend('hto','hte')
%  
%  %4
% figure;
%  semilogy(rd,abs(X(:,4)),'b',rd,abs(Y(:,4)),'r','LineWidth',2); 
%  title('hf') 
%  xlabel('r') 
%  ylabel('|hf|') 
%  legend('hfo','hfe')
 
 
% Графики-----------------------------------------------------------
%1 
 plot(rd,abs(X(:,1)),'b',...
      rd,abs(Y(:,1)),'r','LineWidth',2), grid;
 legend ('et, n0','et, nx')
%2 
figure; plot(rd,abs(X(:,2)),'b',...
        rd,abs(Y(:,2)),'r','LineWidth',2), grid;
 legend ('ef, n0','ef, nx')
%3
figure; plot(rd,abs(X(:,3)),'b',...
        rd,abs(Y(:,3)),'r','LineWidth',2), grid;
legend ('ht, n0','ht, nx')
%4
figure; plot(rd,abs(X(:,4)),'b',...
        rd,abs(Y(:,4)),'r','LineWidth',2), grid;
legend ('hf, n0','hf, nx')

Решение системы ОДУ


Проверяй

*картинку делал так:
Matlab M
>> for i=1:4; subplot(2,2,i),plotyy(rd,abs(X(:,i)),rd,abs(Y(:,i))),grid, legend('x','y'),end
0
0 / 0 / 0
Регистрация: 11.02.2014
Сообщений: 9
12.02.2014, 21:53  [ТС] 3
Графики остались такими же как были, но программа стала приятнее на вид. спасибо!
0
5142 / 3480 / 356
Регистрация: 02.04.2012
Сообщений: 6,387
Записей в блоге: 16
12.02.2014, 23:13 4
А в формулах нет ошибки?
0
0 / 0 / 0
Регистрация: 11.02.2014
Сообщений: 9
14.02.2014, 11:44  [ТС] 5
вроже бы нету ошибок...

Добавлено через 14 часов 31 минуту
Заметил странность: если выводить компоненты тензора Т во время счета программы, то некоторые из них получаются векторами
0
5142 / 3480 / 356
Регистрация: 02.04.2012
Сообщений: 6,387
Записей в блоге: 16
14.02.2014, 18:57 6
Цитата Сообщение от Olegmq Посмотреть сообщение
если выводить компоненты тензора Т во время счета программы, то некоторые из них получаются векторами
просто матрица Т формируется постепенно, поэтому окончательный ее вид получается после присвоения последнего элемента. Если бы некоторые элемнты были бы векторами, то возникала бы ошибка несоответствия размерностей.
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
14.02.2014, 18:57

Заказываю контрольные, курсовые, дипломные работы и диссертации здесь.

Решение системы ОДУ (модель Ресслера)
Доброго времени суток! Не получается реализовать код в матлабе(Сама задача на фотографии). Могли...

Решение неоднородной линейной системы ОДУ в MatLab
Добрый вечер. Необходимо решить неоднородную линейную систему ОДУ. Неоднородность зависит от...

Решение системы ОДУ с условиями на конце интервала
Доброй ночи, форумчани! Подскажите в каком направлении копать. Необходимо решить систему ОДУ: ...

Решение системы ОДУ (стабилизация цены в денежном выражении)
Здравствуйте, Помогите пожалуйста с заданием: x'=a1-a2*x*y, y'=b1-b2*y-b3*x*y При a1=0.3, ...

Решение системы оду, построение графиков и оптимизация кода
Есть система, представленная в м-файле.Грубо говоря,а1..а8 меняются от концентрации лекарства,т.е.:...

Решение ОДУ
Доброго времени суток, имеется задача: x0.5d2y/dx2=y1.5 y(0)=1, y(inf)=0. Попытался решить...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
6
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.