How can I run efficiently a function that changes values in an other function

1 回表示 (過去 30 日間)
Hello everybody,
I am solving an ode with ode45. Here is the equation: xddot + b*xdot + k*x = 9.81. I need to solve it for several values of the pair (k,b) (say more than 10000 pairs)
I wrote this matlab function to be used in the ode45 call:
function Ydot = func(t,Y)
k=2;
b=5;
Ydot=zeros(2,1);
Ydot(1)=Y(2);
Ydot(2)=9.81-b*Y(2)-k*Y(1);
For instance if I call "T=0:1e-3:.5; [t,Y]= ode45('func',T,[0 1]);", it will give in Y, solution (x and xdot) for the equation with (k,b) = (2,5).
(I know that analytic solution can be readily found for this differential equation but this has been just used to illustrate my issue, my real problem is more complex.)
I have a another function that can change the values of k and b to the one provided as arguments. The call of this function is:
NewValuesfunc('func.m',valueOFk,valueOFb);
This function has been tested and it works perfectly. Now I wrote another function, whose inputs are k and b. The function first runs NewValuesfunc to change the values in func.m, second run ode45 with (the new) func.m. The output is only x (not x and xdot). Here is that function:
function y = ComputedResult(x)
T=0:1e-3:.5;
NewValuesfunc('func.m',x(1),x(2)); % change the values first...
[t,Y]= ode45('func',T,[0 1]); %... then find the solution for the new values
y = Y(:,1); % only take x
Now I am testing ComputedResult in a for loop. Initially the values in the func.m file are k=2 and b=5. I save the solution from this pair in xExact (time is always the same: T=0:1e-3:.5). Then I have a set of 4 pairs (k,b) saved in a 2x4 matrix. For each pair ComptedResult is called, which yields the solution (saved in y), this solution is plotted along with the solution from (k,b)=(2,5) (the legend is so that, at each iteration, values of k, b and i are displayed). It turned out that, all the 4 solutions returned by ComputedResult are exactly the solution from (k,b)=(2,5) (initial values), despite the fact that the values in the func.m file has been changed.
Here is the code I used for the test:
T=0:1e-3:.5;
[t,Y]= ode45('func',T,[0 1]); % values are k = 2 and b = 5
xExact=Y(:,1);
kb = [2.5 7.8 8 9;...
.1 2.5 6 4.5]; % a set of 4 pairs of kb, ks in the 1st line and bs in the 2nd
for i=1:4
y = ComputedResult(kb(:,i));
% ComputedResult changes the values of k and b (using kb) in the func.m file first
% and runs ode45 with 'func' to find the corresponding solution of the differential equation
plot(T,xExact,T,y), legend('k=2N/m and b=5N.s/m',sprintf('k =%gN/m and b=%gN.s/m
(i=%d)',kb(1,i),kb(2,i),i)), pause(2);
end
% this should display the curve solution for (k,b)=(2,5) and 4 others different curves (with 2 seconds (pause(2)) of interval between each)
My question is why? Can you tell me please, where the problem is coming from: despite the values in the func.m file has been changing, the computed solution of the differential equation is the same as for the initial pair (k,b)=(2,5).
Thank you in advance for your help.

採用された回答

Walter Roberson
Walter Roberson 2013 年 7 月 4 日
Don't do things this way. Parameterize the call for the ode function. http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

その他の回答 (2 件)

comlich
comlich 2013 年 7 月 4 日
Thank you very much Walter. I just did what is explained on the page you sent me and everything worked very good. Thank you so much!!
However I am still curious about why the way I first did it is not working. So if there is anyone out there who see the reason for that, it will be my pleasure to know it. Thanks

comlich
comlich 2013 年 7 月 4 日
By the way, here is my new ComputedResult function (for those who may be interested):
function y = ComputedResult(x)
T=0:1e-3:.5;
[t,Y]= ode45(@func,T,[0 1]);
function Ydot = func(t,Y)
k=x(1);
b=x(2);
Ydot=zeros(2,1);
Ydot(1)=Y(2);
Ydot(2)=1*9.81-b*Y(2)-k*Y(1);
end
y = Y(:,1);
end

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by