У меня проблема с решением уравнения удержания методом Рунге-Кутта (2-й порядок) в Scilab.
Я должен решить уравнение:
dh/dt=(InF(t)-OutF(t))/F(h), или это уравнение в картинке (ссылка)
where:
h - height [m];
t - time [sec];
InF- inflow to reservoir [m^3/sec];
OutF - outflow from reservoir [m^3/sec];
F - area [m^2];
но у меня нет OutF (t), у меня есть OutF (h), и мне нужно получить приток-отток волны преобразования диаграммы.
Я не знаю, как я могу поместить код scilab без ошибок.
Пример кода (встреч. Рунге-Кутта 2-го порядка) для функции (x, y):
function [y,x]=ExplicitRungeKutta(y0,x0,f,dt,N)
y=zeros(N,1)*%nan;
x=zeros(N,1)*%nan;
y(1)=y0;
x(1)=x0;
for i=2:1:N do
k1=f(x(i-1),y(i-1));
k2=f(x(i-1)+3*dt./4,y(i-1)+3*k1*dt./4);
y(i)=y(i-1)+(k1./3+2*k2./3)*dt;
x(i)=x(i-1)+dt;
end
endfunction
function dy=F(x,y)
dy=-x./(6*y-2*(x.^2)./y);
endfunction
y0=1;
x0=1;
xN=1.34;
N=25;
dt=(xN-x0)/N;
[y,x]=ExplicitRungeKutta(y0,x0,F,dt,N);
plot2d(x,y)
и это хорошая программа, но мой случай сложнее отразить в scilab.
Это код ошибки программы:
clear
clc
WYS=[0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3 3.25 3.5 3.75 4 4.25 4.5 4.75 5 5.25 5.5 5.75 6];
FZ=[7 12.02 22.00 56.03 170.07 346.36 492.35 655.67 820.95 969.25 1164.35 1461.50 1885.19 2429.21 3040.32 3617.32 4155.82 4605.32 5018.49 5448.25 5801.89 6140.23 6414.84 6710.24 7013.57];
//wys - height in metres, Fz - area in m^2
TIME=[0 900 1800 2700 3600 4500 5400 6300 7200 8100 9000 9900 10800 11700 12600 13500 14400 15300 16200 17100 18000 18900 19800 20700 21600 22500 23400 24300 25200 26100 27000 27900 28800 29700 30600 31500 32400 33300 34200 35100 36000 36900 37800 38700 39600 40500 41400 42300 43200 44100 45000 45900 46800 47700 48600 49500 50400 51300 52200 53100 54000 54900 55800 56700 57600 58500 59400 60300 61200 62100 63000 63900 64800 65700 66600 67500 68400 69300 70200 71100 72000 72900 73800 74700 75600 76500 77400 78300 79200 80100 81000 81900 82800 83700 84600 85500 86400 87300 88200 89100 90000 90900 91800 92700 93600 94500 95400 96300 97200 98100 99000 99900 100800 101700 102600 103500 104400 105300 106200 107100 108000 108900 109800 110700 111600 112500 113400 114300 115200 116100 117000 117900 118800 119700 120600 121500 122400 123300 124200 125100 126000 126900 127800 128700 129600 130500 131400 132300 133200 134100 135000 135900 136800 137700 138600 139500 140400 141300 142200 143100 144000 144900 145800 146700 147600 148500 149400 150300 151200 152100 153000 153900 154800 155700 156600 157500 158400 159300 160200 161100 162000 162900 163800 164700 165600 166500 167400 168300 169200 170100 171000 171900 172800]
R10=[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.1 0.1 0.2 0.3 0.3 0.5 0.8 1.3 1.9 2.6 3.3 4.1 4.9 5.7 6.4 7.1 7.7 8.4 9.0 9.5 10.1 10.6 11.0 11.5 11.9 11.7 11.1 10.2 9.0 7.8 6.7 5.8 5.2 4.7 4.4 4.2 4.0 3.9 3.9 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.8 3.9 3.9 3.9 3.9 3.9 3.9 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.1 4.1 4.1 4.1 4.1 4.1 4.1 4.2 3.9 3.6 3.1 2.6 2.0 1.5 1.1 0.8 0.6 0.4 0.3 0.2 0.1 0.1 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
// n=193 //intlow to reservoir, time in seconds, R10- rain 10years, chart inflow to reservoir
WYS1=[0 0.35 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 3.2 3.3 3.4 3.5 4.3 4.5 4.6 4.7 4.8 4.9 5 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6]
oF1=[0 0 0.05 0.15 0.4 0.66 1.1 1.5 2 2.2 2.42 3 3.4 3.8 4.2 4.52 4.84 5.12 5.42 5.65 5.9 6.14 6.37 6.62 6.84 7.04 7.24 7.45 7.65 7.84 8 9.4 9.75 11.11 13.452 16.435 19.93 23.875 28.2 32.915 37.96 41.955 45.95 49.945 53.94 57.935 61.93 65.925]
//outlow from reservoir, of1 - outflow in m3/sec
//y=h,x=t, inf(t=x)=u outf(h=y)=w, Fz(h=y)=z
function [fz,w,u,h,t]=ExplicitRungeKutta(fz0,w0,u0,h0,t0,f,dt,N)
h=zeros(N,1)*%nan;
t=zeros(N,1)*%nan;
h(1)=h0;
t(1)=t0;
u(1)=u0;
w(1)=w0;
fz(1)=fz0;
for i=2:1:N do
k1=f(t(i-1),h(i-1));
k2=f(t(i-1)+3*dt./4,h(i-1)+3*k1*dt./4);
h(i)=h(i-1)+(k1./3+2*k2./3)*dt;
t(i)=t(i-1)+dt;
u(i)=interp1(TIME,R10,t(i),'linear');
w(i)=interp1(WYS,OF1,h(i),'linear');
fz(i)=interp1(WYS,FZ,h(i),'linear');
end
endfunction
function dh=F(u,w,fz)
dh=(u-w)./fz;
endfunction
h0=0;
t0=0;
u0=0;
w0=0;
fz0=7;
tN=172800;
N=17280;
dt=(tN-t0)/N;
[u,w,fz,h,t]=ExplicitRungeKutta(u0,w0,fz0,h0,t0,F,dt,N);
Может быть, у кого-нибудь есть идеи, как это решить?