Initial Condition of Matrix for ODE45 function
Mostrar comentarios más antiguos
Hi All,
I would like to solve simultaneously 2 ode equation, which actually consist of sets of equation. Is there anybody able to explain why the transferred initial condition mw0 couldnt transferred as 2x35 matrix into mw inside the grinding function.
Script as below
clear
%global wsimulmmx
param=[1.1340 0.0336 1.8436 1.4229] ;
%param=fminsearch(@mmxgrind,param0);
psi_lau=param(1)
alfa_lau=param(2)
beta_lau=param(3)
epsilon_lau=param(4)
%wsimulmmx
%function sse=mmxgrind(param);
param;
w0=[0.106 0.126 0.145 0.162 0.18 0.21 0.223 0.249 0.272 0.303 0.355 0.444 0.567 0.731 1.012 1.435 1.897 2.454 3.13 3.683 4.344 5.074 5.825 6.43 7.34 7.584 7.586 7.756 8.031 7.753 6.795 4.206 2.085 1.003 0.504 ];
m0=w0;
mw0=[w0;m0]
[~,w1]=ode45(@grinding,[0 60],mw0,[],param);
[~,w2]=ode45(@grinding,[0 120],mw0,[],param);
[~,w3]=ode45(@grinding,[0 180],mw0,[],param);
[~,w4]=ode45(@grinding,[0 240],mw0,[],param);
[~,w5]=ode45(@grinding,[0 300],mw0,[],param);
[~,w6]=ode45(@grinding,[0 360],mw0,[],param);
[~,w7]=ode45(@grinding,[0 420],mw0,[],param);
[~,w8]=ode45(@grinding,[0 480],mw0,[],param);
wdata= [0 0 0 0 0 0 0.12 0.162 0.183 0.239 0.325 0.452 0.583 0.632 0.873 1.412 1.857 2.403 3.065 3.606 4.253 4.968 5.703 6.266 7.05 7.931 8.029 8.42 8.833 8.561 6.659 3.999 1.902 1.02 0.494 ; % weight fraction pada t=60
0 0 0 0 0 0 0 0 0 0.098 0.181 0.227 0.436 0.553 0.781 1.39 1.828 2.365 3.017 3.549 4.186 4.889 5.613 6.167 6.938 7.806 8.169 8.562 9.448 9.35 7.007 3.955 1.922 1.042 0.521 ; % weight fraction pada t=120
0 0 0 0 0 0 0 0 0 0 0 0.109 0.255 0.408 0.713 1.019 1.427 1.946 2.451 3.159 3.764 4.451 5.198 6.217 7.134 7.950 9.173 9.580 10.386 10.396 6.615 3.950 2.143 1.009 0.547;
0 0 0 0 0 0 0 0 0 0 0 0.109 0.255 0.408 0.713 1.019 1.427 1.946 2.451 3.159 3.764 4.451 5.198 6.217 7.134 7.950 9.173 9.580 10.386 10.396 6.615 3.950 2.143 1.009 0.547;
0 0 0 0 0 0 0 0 0 0 0 0.109 0.255 0.408 0.713 1.019 1.427 1.946 2.451 3.159 3.764 4.451 5.198 6.217 7.134 7.950 9.173 9.580 10.386 10.396 6.615 3.950 2.143 1.009 0.547;
0 0 0 0 0 0 0 0 0 0 0 0.109 0.255 0.408 0.713 1.019 1.427 1.946 2.451 3.159 3.764 4.451 5.198 6.217 7.134 7.950 9.173 9.580 10.386 10.396 6.615 3.950 2.143 1.009 0.547;
0 0 0 0 0 0 0 0 0 0 0 0.109 0.255 0.408 0.713 1.019 1.427 1.946 2.451 3.159 3.764 4.451 5.198 6.217 7.134 7.950 9.173 9.580 10.386 10.396 6.615 3.950 2.143 1.009 0.547;
0 0 0 0 0 0 0 0 0 0 0 0.109 0.255 0.408 0.713 1.019 1.427 1.946 2.451 3.159 3.764 4.451 5.198 6.217 7.134 7.950 9.173 9.580 10.386 10.396 6.615 3.950 2.143 1.009 0.547]; % weight fraction pada t=8 jam
%global wsimul
wsimul= [w1(end,:);
w2(end,:);
w3(end,:);
w4(end,:);
w5(end,:);
w6(end,:);
w7(end,:);
w8(end,:)];
sse=sum(sum(((wdata-wsimul)).^2)')/(35*3) % sums of square error
%end
function dmwdt=grinding(~,mw,param);
%param1=param
%global param %w0 w
F=10 %Flow rate sirkulasi
Vbowl=1000 %Volume bowl
Vmmx=20 %Volume chamber mmx
mw
w=mw(1,:);
m=mw(2,:);
x=[3.905 3.409 2.976 2.599 2.269 1.981 1.729 1.51 1.318 1.151 1.005 0.877 0.766 0.669 0.584 0.51 0.445 0.389 0.339 0.296 0.259 0.226 0.197 0.172 0.15 0.131 0.115 0.1 0.087 0.076 0.067 0.058 0.051 0.044 0.039 ];
for k=1:length(x)
S(k)=(param(2)*x(k)^param(3))/(1+x(k)^param(4));
for u=1:k
if k==1
deltaB(k,u)=((x(k)*1.145)/x(u))^param(1)-((x(k)/x(u))^param(1));
else
deltaB(k,u)=((x(k-1))/x(u))^param(1)-((x(k)/x(u))^param(1));
end
end
end
for k=1:length(x);
WStemp(k)=m(k).*S(k);
end
WS=WStemp;
deltaB;
summation=WS*deltaB';
%pause
for k=1:length(x);
dmdt(k)=F/Vmmx*(w(k)-m(k))+summation(k)-S(k)*m(k);
dwdt(k)=F/Vbowl*(m(k)-w(k));
end
dmdt=dmdt';
dwdt=dwdt';
dmwdt=[dwdt;dmdt];
end
Result as below
>> mmx
psi_lau =
1.1340
alfa_lau =
0.0336
beta_lau =
1.8436
epsilon_lau =
1.4229
mw0 =
Columns 1 through 13
0.1060 0.1260 0.1450 0.1620 0.1800 0.2100 0.2230 0.2490 0.2720 0.3030 0.3550 0.4440 0.5670
0.1060 0.1260 0.1450 0.1620 0.1800 0.2100 0.2230 0.2490 0.2720 0.3030 0.3550 0.4440 0.5670
Columns 14 through 26
0.7310 1.0120 1.4350 1.8970 2.4540 3.1300 3.6830 4.3440 5.0740 5.8250 6.4300 7.3400 7.5840
0.7310 1.0120 1.4350 1.8970 2.4540 3.1300 3.6830 4.3440 5.0740 5.8250 6.4300 7.3400 7.5840
Columns 27 through 35
7.5860 7.7560 8.0310 7.7530 6.7950 4.2060 2.0850 1.0030 0.5040
7.5860 7.7560 8.0310 7.7530 6.7950 4.2060 2.0850 1.0030 0.5040
F =
10
Vbowl =
1000
Vmmx =
20
mw =
0.1060
0.1060
0.1260
0.1260
0.1450
0.1450
0.1620
0.1620
0.1800
0.1800
0.2100
0.2100
0.2230
0.2230
0.2490
0.2490
0.2720
0.2720
0.3030
0.3030
0.3550
0.3550
0.4440
0.4440
0.5670
0.5670
0.7310
0.7310
1.0120
1.0120
1.4350
1.4350
1.8970
1.8970
2.4540
2.4540
3.1300
3.1300
3.6830
3.6830
4.3440
4.3440
5.0740
5.0740
5.8250
5.8250
6.4300
6.4300
7.3400
7.3400
7.5840
7.5840
7.5860
7.5860
7.7560
7.7560
8.0310
8.0310
7.7530
7.7530
6.7950
6.7950
4.2060
4.2060
2.0850
2.0850
1.0030
1.0030
0.5040
0.5040
Index exceeds matrix dimensions.
Sorry for long script,
Thanks
Respuestas (1)
darova
el 16 de Mzo. de 2019
Try this
function dmwdt=grinding(t,mw,param);
%% ...
ode45(@(t,mw)@grinding(t,mw,param),[0 60],mw0);
3 comentarios
Khanif Eko Prasetyo
el 16 de Mzo. de 2019
darova
el 16 de Mzo. de 2019
You changed your code according my suggestion and it triggers error? Can you show an exact error?
Khanif Eko Prasetyo
el 22 de Mzo. de 2019
Categorías
Más información sobre Programming en Centro de ayuda y File Exchange.
el 16 de Mzo. de 2019
el 22 de Mzo. de 2019
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
