Changing parameters in an ODE
4 次查看(过去 30 天)
显示 更早的评论
I wish to change one of the parameters at each time step for an ODE solution plot. My function is as follows:
function [t,v]=shig(b,p,m,yo)
[t v] = ode45(@fnsirtry,[0 12],yo);
function fnsir = fnsirtry(t,v)
a = 0.25;
r = 0.14;
fnsir(1) = p - m*v(1)-b*v(1)*v(2)+a*v(3);
fnsir(2) = (b*v(1))- (m + r)*v(2);
fnsir(3) = (r*v(2))-((m+a)*v(3));
fnsir = fnsir(:);
end
end
When I plot the ODE, I want to change the value of b at each step. I have used a code like teh following:
temp = [5 9 12 17 19 24 28 27 22 17 10 7];
beta = 0.0000025*temp;
i = [1 2 3 4 5 6 7 8 9 10 11 12];
for i = 1;
b = beta(:,1)
end
for i = 2;
b = beta(:,2)
end
for i = 3;
b = beta(:,3)
end
for i = 4;
b = beta(:,4)
end
for i = 5;
b = beta(:,5)
end
for i = 6;
b = beta(:,6)
end
for i = 7;
b = beta(:,7)
end
for i = 8;
b = beta(:,8)
end
for i = 9;
b = beta(:,9)
end
for i = 10;
b = beta(:,10)
end
for i = 11;
b = beta(:,11)
end
for i = 12;
b = beta(:,12)
end
P = 5000;
m = 0.013;
yo = [200000 160 0];
[t v] = shig(b,p,m,yo);
w = [159
148
143
137
102
91
85
137
119
108
104
100];
subplot (1,2,1)
plot(t,v(:,2))
%,'-r*','Linewidth',1.5,'MarkerSize',5)
title('Infected Population')
%legend('Disease free state','Test state')
subplot(1,2,2)
plot(w)
Could I please get some light on this,, pls...
采纳的回答
Star Strider
2015-8-30
I don’t understand all the for loops. Assuming your ODE works and integrates as you want it to (I didn’t run your code), I would just do:
for k1 = 1:12
[t{k1}, v{k1}] = shig(beta(:,k1),p,m,yo);
end
and then plot the individual cell vectors. Note that if you define your evaluation times as a vector of discrete times for all integrations, rather than as a two-element range, you can use a matrix to store them rather than a cell array. Your call.
20 个评论
Ojaswita
2015-8-30
Thanks alot for the help. Yes my code runs and integrates appropriately. Could you please explain how to plot individual cell vectors and how to define the evaluation times as a vector.
You have understood what I wanted to do, thanks alot! I tried using your piece of code, but I get the following error:
??? Comma separated list expansion has cell syntax for an array that is not a cell.
Star Strider
2015-8-30
My pleasure.
First, to use ‘tspan’ (as referred to in the ODE solver documentation) as a vector, just define it as such:
tspan = linspace(0, 12, 50);
creates a 50-element vector going from 0 to 12.
Second, I don’t know precisely what the problem is, but see if changing the loop to this version solves it:
for k1 = 1:12
[tv, vv] = shig(beta(:,k1),p,m,yo);
t{k1} = tv;
v{k1} = vv;
end
The sort of assignment I used initially doesn’t usually cause problems.
Third, to plot individual cell vectors, considering you seem to want to plot them all as subplots, might be something like this:
for k1 = 1:12
subplot(12, 1, k1)
plot(t{k1}, v{k1}(:,2))
grid
end
This is untested code as well, so I leave you to experiment with it to see if it works.
Ojaswita
2015-8-30
编辑:Ojaswita
2015-8-30
Thanks for the help... the errors still persist on the for loop... I also cant figure out what exactly the problem is as it looks fine to me. I would greatly appreciate if you could try to run my code... I have ti get a separate value of b for each step...
I tried to use parenthesis instead of curly brackets, and I got the following error:
??? In an assignment A(I) = B, the number of elements in B and I must be the same.
Star Strider
2015-8-30
The parentheses will work if you create a fixed-length time vector:
b = [3 13]; % Proxy For ‘beta’ Array
tv = linspace(0, 2, 75); % Define Fixed-Length Evaluation Times (NECESSARY)
for k1 = 1:2
ode = @(t,x) [x(1); x(1) - x(2)*sin(b(k1).*pi.*t/2)]; % Define ODE Function With Changing ‘b’-Element
[t(:,k1),v(:,:,k1)] = ode45(ode, tv, [1; 1]); % Integrate ODE
end
figure(1)
for k1 = 1:length(b)
subplot(2,1,k1)
plot(t(:,k1), v(:,2,k1))
grid
end
I created a simple ODE for this, but yours should work with it. Be sure to define tspan as a vector and not a range, or this will fail.
I don’t have your constants so I can’t run your code. However this should work with your code, since it does in my simulation. The cell arrays also work in my simulation, so I have no idea why they don’t in your code.
Ojaswita
2015-8-31
编辑:Ojaswita
2015-8-31
Im getting confused badly... Im trying to put my work in the above format and getting errors... Could you try my code, that would be extremely extremely helpful. My constants are: p = 5000; m =0.013; a =0.25; r = 0.14;
My code (according to the above format) is:
temp = [5 9 12 17 19 24 28 27 22 17 10 7];
b = 0.0000025*temp;
tv = linspace(0,12,50);
yo = [200000 160 0];
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;
for k1 = 1:12
ode = @(t,v) [p - m*v(1)-b*v(1)*v(2)+a*v(3); (b*v(1))- (m + r)*v(2);(r*v(2))-((m+a)*v(3));];
[t(:,k1),x(:,:,k1)] = ode45(ode, tv, yo);
end
figure(1)
for k1 = 1:length(b)
subplot(2,1,k1)
plot(t(:,k1), x(:,2,k1))
grid
end
-----------------------------------------------------
And when i tried to run the previous assignment format, like this:
temp = [5 9 12 17 19 24 28 27 22 17 10 7];
beta = 0.0000025*temp;
P = 5000;
m = 0.013;
yo = [200000 160 0];
for k1 = 1:12
b = beta(:,k1);
[t{k1}, v{k1}] = shig(b,p,m,yo);
end
subplot (1,2,1)
plot(t,v(:,2))
I get an error that says:
??? Error using ==> plot Size mismatch in Param Cell / Value Cell pair Error in ==> compare at 64 plot(t,v(:,2))
Star Strider
2015-8-31
Since you have 12 values of ‘b’, you have to run the integration loop 12 times, each with a different value for ‘b’. Then create 12 subplots.
This works:
temp = [5 9 12 17 19 24 28 27 22 17 10 7];
b = 0.0000025*temp;
tv = linspace(0,12,50);
yo = [200000 160 0];
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;
for k1 = 1:12
ode = @(t,v) [p - m*v(1)-b(k1)*v(1)*v(2)+a*v(3); (b(k1)*v(1))- (m + r)*v(2);(r*v(2))-((m+a)*v(3));];
[t(:,k1),x(:,:,k1)] = ode45(ode, tv, yo);
end
NrPlots = length(b);
figure(1)
for k1 = 1:NrPlots
subplot(NrPlots,1,k1)
plot(t(:,k1), x(:,2,k1))
grid
end
Note that in this code, I subscripted ‘b’ as b(k1) in your ODE. That, and adjusting the subplots accordingly were the only changes needed.
Ojaswita
2015-8-31
First of all, I really appreciate your assistance. Thank you so very much!! Its really encouraging!
I do not need subplots... I think we lost track of each other there...
I want to plot the second graph of the solution of my ODE, for twelve points, with each point having a different value of b in the solution.
Star Strider
2015-8-31
My pleasure.
I’m not quite sure what you mean by ‘twelve points, with each point having a different value of b in the solution.’ The code produces twelve solutions.
See if this does what you want:
figure(2)
plot(t, squeeze(x(:,2,:)))
lgndstr = regexp(sprintf('\\beta = %.3E\n',b), '\n', 'split');
legend(lgndstr(1:end-1), 'Location','SW', 'FontSize',6)
grid
Ojaswita
2015-9-1
Ok, ill try to be more accurate. I want to change the parameter value of b at every time interval in the ODE solution.
I want to plot the second equation of the ODE solution for 12 time intervals. And at each time interval, I want b to have a different value for the solution.
I hope im more clear...
Star Strider
2015-9-1
I’m lost. Do you want the ODE solved at 12 points only, changing ‘b’ at each time point? (I thought that’s what I did.)
How do you define the time in that instance? How do you define the times when ‘b’ changes?
If you want to define a time interval and then change ‘b’ at specific times, that’s relatively straightforward. I just need specifically to know what you want.
Ojaswita
2015-9-1
编辑:Ojaswita
2015-9-1
I tried to further ammend my code as per my understanding but now it appears to be doing the initial thing using only the first value of b.
Here is the function:
function [t,v]=shig(b,p,m,yo)
[t v] = ode45(@fnsirtry,tv,yo);
function fnsir = fnsirtry(t,v)
temp = [5 9 12 17 19 24 28 27 22 17 10 7];
beta = 0.0000025*temp;
a = 0.25;
r = 0.14;
for k=1:12
b = beta(:,k)
fnsir(1) = p - m*v(1)-b*v(1)*v(2)+a*v(3);
fnsir(2) = (b*v(1))- (m + r)*v(2);
fnsir(3) = (r*v(2))-((m+a)*v(3));
fnsir = fnsir(:);
end
end
end
And the plotting:
t = [5 9 12 17 19 24 28 27 22 17 10 7];
b = 0.0000025*t
tv = linspace(0,12,50);
yo = [200000 160 0];
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;
for k = 1:12
b(k) = B(:,k)
[t v]=shig(b(k),p,m,yo)
end
subplot(1,2,1)
plot(t,v(:,2))
Now it appears to be taking only the first value...
Star Strider
2015-9-1
Please answer the questions I asked in my previous Comment. I need to understand what you’re doing.
Ojaswita
2015-9-1
Sorry, I missed one of your comments.
Yes I need the ODE solved at 12 points only, changing ‘b’ at each time point. But I do not want 12 separate subplots, that is what I got. The initial assignment technique made sense to me also, but its giving an error...
b is dependent on temperature, which i have defined as temp. temp gives the average temperature per month. b is a constant value multiplied by the matrix temp.
I need to plot this graph to compare it with original data so that I manage to fix b to give the most perfect match.
Star Strider
2015-9-1
See if this does what you want:
temp = [5 9 12 17 19 24 28 27 22 17 10 7];
b = 0.0000025*temp;
% tv = linspace(0,12,50);
time = 1:12;
yo = [200000 160 0];
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;
for k1 = time
tv = [k1-0.5, k1, k1+0.5]; % Three-Element Time Vector
ode = @(t,v) [p - m*v(1)-b(k1)*v(1)*v(2)+a*v(3); (b(k1)*v(1))- (m + r)*v(2);(r*v(2))-((m+a)*v(3));];
[t(:,k1),v(:,:,k1)] = ode45(ode, tv, yo);
yo = v(3,:,k1); % Updated Initial Conditions Vector For Next Iteration
vp(k1) = yo(2); % Element To Be Plotted
end
figure(2)
plot(time(1), vp(1), 'p')
hold on
for k1 = 2:length(time)
plot(time(k1), vp(k1), 'p')
end
hold off
lgndstr = regexp(sprintf('\\beta = %.3E\n',b), '\n', 'split');
legend(lgndstr(1:end-1), 'Location','SW', 'FontSize',6)
grid
I created a three-element time vector for each integration, and changed ‘b’ with each. The previous end conditions became the new initial conditions. This is the standard way of ‘piecewise-integrating’ an otherwise discontinuous differential equation system. I added ‘vp’ so you could specify it correctly if I got it wrong. It still saves all the results in the ‘v’ matrix.
Ojaswita
2015-9-3
THANK YOU!!! At least its absolutely somewhere what I needed!
But its still abit different... And im not getting confused in my own mesh... Can we have a conversation out of the forum, please please...
Star Strider
2015-9-3
My pleasure!
I guess I’m still not understanding exactly what you want to do. What about it is ‘a bit different’ from what you need? How is it different?
I prefer to keep our conversations here because Answers allows markup to make code easier to read, and that may be important. I also can refer to everything else I wrote here, and I know where to look to find all posts in this thread.
Ojaswita
2015-9-7
I managed to crack my head and get the issue resolved - thank you so much. I was doing a mistake by adding the b definition in both my files hence the code was not working.
But I need some further help with this as the function is complicated. I used a one matrix and an easy formula just to get the code working. Here is the code i used:
function [t,v]=shig(b,p,m,yo)
tv = linspace(0,11,5);
[t v] = ode45(@fnsirtry,tv,yo);
function fnsir = fnsirtry(t,v)
a = 0.25;
r = 0.14;
fnsir(1) = p - m*v(1)-b*v(1)*v(2)+a*v(3);
fnsir(2) = (b*v(1))- (m + r)*v(2);
fnsir(3) = (r*v(2))-((m+a)*v(3));
fnsir = fnsir(:);
end
end
and...
temp = [5 9 12 17 19 24 28 27 22 17 10 7];
beta = 0.0025*temp;
tv = linspace(0,12,5);
yo = [200000 160 0];
p = 5000;
m = 0.013;
a = 0.25;
r = 0.14;
for k = 1:12
b = beta(k)
[t,v]=shig(b,p,m,yo)
end
subplot(1,2,1)
plot(t,v(:,2))
But I need amendments.. My b is actually a piecewise function.
whereby for 0-6, it is 120*cos(k) for 6-7, it is 50*k - 265 for 7-11 it is 10^k +7
so I need to adjust that.. .would appreciate if you could further help me... And I thank you so much for your help.
Star Strider
2015-9-7
If it’s piecewise, you have to integrate it over each piece with a different call to ode45.
Do not integrate over the discontinuities! No ode solver does that well.
Use the previous final values of your variables as the initial conditions for the next values of them (or choose a different set, if that works for you). I used the previous values in the code in my comment of 1 Sep 2015 at 17:19, but that is essentially the way to do what you want. It depends on how discontinuous your function is. Note that you have to update your independent variable as well, so it goes from 0-6, 6-7, 7-11. Be sure the end value of the previous interval doesn’t overlap with the first value of the next interval.
Ojaswita
2015-9-10
Yes, I get it. Thanks so much for the help! I'll do it and come back if I face troubles. Thanks a ton for all the time help! I highly appreciate it. Students like me look forward to such sites for help so that we may learn more and apply better. Thank you!
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
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 (한국어)