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

Answer 1

[darova]

Try this

function dmwdt=grinding(t,mw,param);
%% ...
ode45(@(t,mw)@grinding(t,mw,param),[0 60],mw0);

[Khanif Eko Prasetyo]

Unfortunately it doesnt work, it trigger parenthesis error

[darova]

You changed your code according my suggestion and it triggers error? Can you show an exact error?

[Khanif Eko Prasetyo]

eventually I do the old fashion way to solve the issue, by rearranging the mw into required matrix dimension manually. Thanks anyway for your answer