Unable to perform assignment because the left and right sides have a different number of elements.
1 次查看(过去 30 天)
显示 更早的评论
Hi all
When I run the attached code, I get the error message:
Unable to perform assignment because the left and right sides have a different number of elements.
Error in DynamicOptimization10Jan2019>myModel (line 90)
dxdt(1)= -2.85.*(1/((1/0.00639039)+(1/k)).*12.54.*x.*x.*(1-(2.69636.*x)/(1+(2.69636.*x))));
Error in DynamicOptimization10Jan2019>@(t,x)myModel(t,x,k) (line 81)
[t,x] = ode15s(@(t,x)myModel(t,x,k),time,x0);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 150)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in DynamicOptimization10Jan2019>myObj (line 81)
[t,x] = ode15s(@(t,x)myModel(t,x,k),time,x0);
Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});
Error in DynamicOptimization10Jan2019 (line 3)
k = fminsearch(@myObj,k0);
I know the error I am getting implies that I am trying to put more elements into less elements, or vice versa. However, I do not know how to make elements equal.
please help.
Kind Regards
16 个评论
KSSV
2019-1-10
YOur output in the line:
dxdt(1)= -2.85.*(1/((1/0.00639039)+(1/k)).*12.54.*x.*x.*(1-(2.69636.*x)/(1+(2.69636.*x))));
is a matrix..and you are trying to save it as array. You need to rethink on your code.
Dursman Mchabe
2019-1-10
Thanks a lot such a quick response. I will rethink. My first move will be to replace .* with *.
Dursman Mchabe
2019-1-10
Oops, it gives another error:
Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To perform
elementwise multiplication, use '.*'.
Error in DynamicOptimization10Jan2019>myModel (line 89)
dxdt(1)= -2.85*(1/((1/0.00639039)+(1/k))*12.54*x*x*(1-(2.69636*x)/(1+(2.69636*x))));
I must rethink.
Dursman Mchabe
2019-1-10
Thanks a lot for the comment Walter. The challenge is that I have 4 ODEs dxdt(1) ... dxdt(4), with only 2 varibles, x(1) and x(2). Perhaps I should try to pass x as x(1) and x(2).
Dursman Mchabe
2019-1-10
Oops, it gives:
Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-4.
Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});
Error in DynamicOptimization10Jan2019 (line 3)
k = fminsearch(@myObj,k0);
How can I change the left to a 1-by-4?
Walter Roberson
2019-1-10
You cannot change the left to 1 x 4. You are using myObj as the objective function for fminsearch, and objective functions for fminsearch must always return scalars. Returning 1 x 4 there would be like trying to simultaneously minimize 4 things, but without instructions about the relative values of decreases between the four parts (because you could easily encounter a situation where a small increase in one of the four is needed to permit a large decrease in one or more of the other three.)
If you do want simultaneous minimization of four outputs, then you should look at gamultiobj and examine the pareto front.
Dursman Mchabe
2019-1-10
I want to try and keep things simple. When I hypothetically solve just one ODE, the code works fine. See below:
%% solve with fminsearch
k0 = 59563518.8216;
k = fminsearch(@myObj,k0);
disp(['fminsearch: k = ' num2str(k)]);
mySim(k);
%% solve with fmincon
LB = 20;
UB = 100;
k = fmincon(@myObj,k0,[],[],[],[],LB,UB);
disp(['fmincon: k = ' num2str(k)]);
mySim(k);
%% solve graphically
n = 100;
k = linspace(4e8,5e9,n);
for i = 1:n,
obj(i) = myObj(k(i));
end
figure(2)
plot(k,obj,'b--','LineWidth',3);
xlabel('k value')
ylabel('Objective Value')
%% confidence interval
%(S(k) - S(K)) / S(K) <= p / (n-p) * F(p,n-p,1-alpha)
p = 1; % number of parameters
np = 30; % number of data points
alpha = 0.05; % alpha, confidence interval
rhs = p / (np-p) * finv(1-alpha,p,(np-p));
ub = min(obj) * rhs + min(obj);
hold on
plot([min(k) max(k)],[ub ub],'r-','LineWidth',3);
legend('Objective','Confidence Limit')
% axis([0.06 0.09 0 1])
% axis([0.07 0.09 0 1])
axis([3e8 5e8 0 7e9])
axis([4e8 5e8 0 7e9])
function obj = myObj(k)
% measured values
time = [1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
];
y = [9.54992586
5.623413252
2.884031503
0.489778819
0.048977882
0.046773514
0.032359366
0.022908677
0.016595869
0.014125375
0.011748976
0.01
0.00851138
0.007585776
0.00691831
0.006309573
];
% predicted values
x0 = 9.54992586;
[t,x] = ode15s(@(t,x)myModel(t,x,k),time,x0);
% compute sum of squared errors
obj = sum((x-y).^2);
end
function dxdt = myModel(t,x,k)
dxdt= -2.85*(1/((1/0.00639039)+(1/k))*12.54*x*x*(1-(2.69636*x)/(1+(2.69636*x))));
end
function [] = mySim(k0)
% Predicted values
x0 = 9.54992586;
[t,x] = ode15s(@(t,x)myModel(t,x,k0),[0 31],x0);
% Actual values
time = [1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
];
y = [9.54992586
5.623413252
2.884031503
0.489778819
0.048977882
0.046773514
0.032359366
0.022908677
0.016595869
0.014125375
0.011748976
0.01
0.00851138
0.007585776
0.00691831
0.006309573
];
figure(1)
hold off
plot(t,x)
hold on;
plot(time,y,'ro')
xlabel('Time')
ylabel('Value')
legend('Predicted','Actual')
end
Is there a better way of applying this code on more than one ODEs?
Walter Roberson
2019-1-10
It is pretty late where I am, but at the moment I get the impression that this is effectively a boundary value problem, trying to find the k that leads to a particular outcome.
I am also suspecting that your model has a closed form solution if approached analytically, so I not certain you need the ode15s call. However, I am too tired to figure out what your original equations must have been.
Dursman Mchabe
2019-1-10
I'm sorry to keep you up. I truly appreciate your help. Where I am it is 11:32 am.
You are absolutely right. I am looking for a k that will make model calculations equal measured values.
Dursman Mchabe
2019-1-10
When I take lot of gueses manually, I find k0 = 46.83370 to work see below:
function [] = mySim(k0)
k0 = 46.83370;
% Predicted values
x0 = [9.54992586;19.89;0;0];
[t,x] = ode15s(@(t,x)myModel(t,x,k0),[0 31],x0);
% Actual values
time = [1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
];
y = [9.54992586 19.89 0 0
5.623413252 18.51227628 1.377723722 0.000176583
2.884031503 17.5510897 2.338910301 0.000584451
0.489778819 16.71100104 3.178998962 0.004671909
0.048977882 16.55633404 3.333665957 0.048352564
0.046773514 16.55556058 3.33443942 0.050608503
0.032359366 16.55050298 3.339497016 0.072770512
0.022908677 16.54718695 3.342813047 0.101976407
0.016595869 16.54497193 3.345028067 0.139244055
0.014125375 16.54410509 3.345894907 0.16245708
0.011748976 16.54327127 3.346728731 0.193464978
0.01 16.54265759 3.347342407 0.225067572
0.00851138 16.54213527 3.34786473 0.26139843
0.007585776 16.5418105 3.348189503 0.290553969
0.00691831 16.5415763 3.348423702 0.315962812
0.006309573 16.54136271 3.348637294 0.34334241];
figure(1)
hold off
plot(t,x)
hold on;
plot(time,y,'v')
xlim([0 35])
ylim([-1 20])
xlabel('Time (sec)')
ylabel('Concentration (mol/m^3)')
legend('Exp-H^+', 'Exp-CaCO_3', 'Exp-Ca^2+', 'Exp-HCO^-_3','Mod-H^+', 'Mod-CaCO_3', 'Mod-Ca^2+', 'Mod-HCO^-_3', 'Location','best')
end
function dxdt = myModel(t,x,k)
dxdt=zeros(4,1);
dxdt(1)= -2.85*(1/((1/0.00639039)+(1/k))*12.54*x(2)*x(1)*(1-(2.69636*x(1))/(1+(2.69636*x(1)))));
dxdt(2)= -(1/((1/0.00639039)+(1/k))*12.54*x(2)*x(1)*(1-(2.69636*x(1))/(1+(2.69636*x(1)))));
dxdt(3)= (1/((1/0.00639039)+(1/k))*12.54*x(2)*x(1)*(1-(2.69636*x(1))/(1+(2.69636*x(1)))));
dxdt(4)= (1/((1/0.00639039)+(1/k))*12.54*x(2)*x(1)*(1-(2.69636*x(1))/(1+(2.69636*x(1)))));
end
How ever I need to use an objective function for a lot other similar codes.
Walter Roberson
2019-1-10
can you post the equation form of the differentiation equationss?
You should also look at bvp4c and related functions .
Dursman Mchabe
2019-1-11
I have found partial success when using lsqcurvefit. It can estimate k. Yet I don't know how to do sensivity analysis and to determine confidence interval without the objective function. Please see the code below.
function KineticModelat30deg
function C=KineticModel(k,t)
c0=[9.54992586;19.89;0;0];
[T,Cv]=ode45(@ODEs,t,c0);
%
function dC=ODEs(t,c)
dcdt=zeros(4,1);
dcdt(1)= -2.85*(1/((1/0.00639039)+(1/k(1)))*12.54*c(2)*c(1)*(1-(k(2)*c(1))/(1+(k(2)*c(1)))));
dcdt(2)= -(1/((1/0.00639039)+(1/k(1)))*12.54*c(2)*c(1)*(1-(k(2)*c(1))/(1+(k(2)*c(1)))));
dcdt(3)= (1/((1/0.00639039)+(1/k(1)))*12.54*c(2)*c(1)*(1-(k(2)*c(1))/(1+(k(2)*c(1)))));
dcdt(4)= (1/((1/0.00639039)+(1/k(1)))*12.54*c(2)*c(1)*(1-(k(2)*c(1))/(1+(k(2)*c(1)))));
dC=dcdt;
end
C=Cv;
end
t=[1
3
5
7
9
11
13
15
17
19
21
23
25
27
29
31
];
c=[9.54992586 19.89 0 0
5.623413252 18.51227628 1.377723722 0.000176583
2.884031503 17.5510897 2.338910301 0.000584451
0.489778819 16.71100104 3.178998962 0.004671909
0.048977882 16.55633404 3.333665957 0.048352564
0.046773514 16.55556058 3.33443942 0.050608503
0.032359366 16.55050298 3.339497016 0.072770512
0.022908677 16.54718695 3.342813047 0.101976407
0.016595869 16.54497193 3.345028067 0.139244055
0.014125375 16.54410509 3.345894907 0.16245708
0.011748976 16.54327127 3.346728731 0.193464978
0.01 16.54265759 3.347342407 0.225067572
0.00851138 16.54213527 3.34786473 0.26139843
0.007585776 16.5418105 3.348189503 0.290553969
0.00691831 16.5415763 3.348423702 0.315962812
0.006309573 16.54136271 3.348637294 0.34334241
];
k0=[0.1;50];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@KineticModel,k0,t,c);
fprintf(1,'\tRate and Adsorption constants:\n')
for i = 1:length(k)
fprintf(1, '\t\tk(%d) = %8.5f\n', i, k(i))
end
tv = linspace(min(t), max(t));
Cfit = KineticModel(k, tv);
figure(1)
plot(t, c, 'v')
hold on
plot(tv, Cfit);
hold off
xlabel('Time (sec)')
ylabel('Concentration (mol/m^3)')
legend('Exp-H^+', 'Exp-CaCO_3', 'Exp-Ca^2+', 'Exp-HCO^-_3','Mod-H^+', 'Mod-CaCO_3', 'Mod-Ca^2+', 'Mod-HCO^-_3', 'Location','E')
end
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Performance and Memory 的更多信息
标签
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)