Running ode45 from problem statement
1 回表示 (過去 30 日間)
古いコメントを表示
ok so my problem statement is "Use the MATLAB built-in function ode45 to perform the same calculations and generate the plots (study Example 10-11)." I have attatched the picture of the problem statement as well and I was able to run the " Sys2ODEsRK4" No problem but Im having trouble with the ode45 :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/285231/image.jpeg)
This is the Sys2ODEsRK4 code
function [t,y,ydot] = Sys2ODEsRK4(ODE1,ODE2,a,b,h,x1,y1)
t(1) = a; y(1) = .025; ydot(1) = 0;
n= (b-a)/h;
for i=1:n
t(i+1) = t(i) + h;
tm = t(i) + h/2;
Kx1 = ODE1(t(i),y(i),ydot(i));
Ky1 = ODE2(t(i),y(i),ydot(i));
Kx2 = ODE1(tm,y(i)+ Kx1*h/2,ydot(i)+Ky1*h/2);
Ky2 = ODE2(tm,y(i)+ Kx1*h/2,ydot(i)+Ky1*h/2);
Kx3 = ODE1(tm,y(i)+ Kx2*h/2,ydot(i)+Ky2*h/2);
Ky3 = ODE2(tm,y(i)+ Kx2*h/2,ydot(i)+Ky2*h/2);
Kx4 = ODE1(t(i + 1),y(i)+ Kx3*h,ydot(i)+Ky3*h);
Ky4 = ODE2(t(i + 1),y(i)+ Kx3*h,ydot(i)+Ky3*h);
y(i+1) = y(i)+(Kx1 + 2*Kx2 + 2*Kx3 + Kx4)*h/6;
ydot(i+1) = ydot(i)+(Ky1 + 2*Ky2 + 2*Ky3 + Ky4)*h/6;
end
Now for my script file part two I got the correct graphs
clc; clear all
a=0;b=10;h=0.02;
[t,y,ydot]=Sys2ODEsRK4(@dydt,@dydotdt,a,b,h,0.025,0);
figure(1)
plot(t,y,'LineWidth',2)
xlabel('Time(s)')
ylabel('Displacement (m)')
title('Displacement over time')
grid on
figure(2)
plot(t,ydot,'LineWidth',2)
xlabel('Time(s)')
ylabel('Velocity (m/s)')
title('Velocity vs Time')
grid on
function dydx=dydt(t,y,ydot)
dydx=ydot;
end
function dydx=dydotdt(t,y,ydot)
L=0.1;g=9.80;
dydx=(-2*g/L)*y;
end
Now when I do my ode45 part this is where im getting my error in the line 6 of this code:
%% Solution using matlab ode45() function
timespan = [0 10];
initCond = 0;
[t,y,ydot] = ode45(@dydt,@dydotdt,timespan,initCond);
%% plotting results
figure
plot(t,y(:,1));
xlabel('time(s)');
ylabel('Displacement (m)');
title('ode45: Displacement vs time');
figure
plot(t,y(:,2));
xlabel('time(s)');
ylabel('Velocity (m/s)');
title('ode45: Velocity vs time');
function dydx=dydt(t,y,ydot)
dydx=ydot;
end
function dydx=dydotdt(t,y,ydot)
L=0.1;g=9.80;
dydx=(-2/L*g)*y;
end
Now this is the example in the textbook for ode45 if this helps but Im lost on what im doing wrong:
function dNdt = PopRate(t,N)
bG = 1.1; bL = 0.00025; dG = 0.005;dL = 0.7;
f1 = bL*N(1)*N(2)-dL*N(1);
f2 = bG*N(2) - dG*N(2)*N(1);
dNdt = [f1;f2];
tspan = [0 25];
Nini = [ 500 3000];
[Time Pop] = ode45(@PopRate, tspan, Nini);
plot(Time,Pop(:,1),'-',Time,Pop(:,2),'--')
xlabel('Time (yr)')
ylabel('Population')
legend('Lions','Gazelles')
0 件のコメント
採用された回答
Peter O
2020 年 4 月 17 日
Hi Daniel,
For your problem, you were very close and had the solutions to each equation worked out correctly. The function signature for ode45 does not spit back the differential. This is a situation where you can take advantage of formatting the differential equations as a vector and have MATLAB solve them simultaneously. Essentially, you have a two state system, where one state is the derivative of the other, but no state needs to be integrated more than once at any timestep. That is:
dx = f(x), where x is a 2-element vector having x(1) equal to y and x(2) as y_dot. Then the solution becomes:
y = x(1);
ydot = x(2);
dx(1) = ydot;
dx(2) = -2gy/L
This is a very powerful trick that is used routinely in solving complex systems.
%% Solution using matlab ode45() function
timespan = [0 10];
initCond = [0.025, 0]; % y, ydot, grouped together
[t,y] = ode45(@dydt,timespan,initCond);
%% plotting results
figure
plot(t,y(:,1));
xlabel('time(s)');
ylabel('Displacement (m)');
title('ode45: Displacement vs time');
figure
plot(t,y(:,2));
xlabel('time(s)');
ylabel('Velocity (m/s)');
title('ode45: Velocity vs time');
function dx=dydt(t,x)
% Compute both derivatives. The states (y and ydot) come in as the x vector, of size 2x1 for time t
L=0.1;
g=9.80;
y = x(1);
ydot = x(2);
dx = zeros(2,1);
dx(1) = ydot;
dx(2) =(-2/L*g)*y;
end
4 件のコメント
Peter O
2020 年 4 月 17 日
It seems to run OK here from a single file using MATLAB's ability to nest a function within a script. Is it potentially calling an old version of dydt? The variable x is undefined in your original code. Type
which dydt
and make sure it's executing the file you expect.
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!