Initial Condition of Matrix for ODE45 function
이전 댓글 표시
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
답변 (1개)
darova
2019년 3월 16일
Try this
function dmwdt=grinding(t,mw,param);
%% ...
ode45(@(t,mw)@grinding(t,mw,param),[0 60],mw0);
댓글 수: 3
Khanif Eko Prasetyo
2019년 3월 16일
darova
2019년 3월 16일
You changed your code according my suggestion and it triggers error? Can you show an exact error?
Khanif Eko Prasetyo
2019년 3월 22일
카테고리
도움말 센터 및 File Exchange에서 Programming에 대해 자세히 알아보기
2019년 3월 16일
2019년 3월 22일
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
