Runge-Kutta 4th Order

11 次查看（过去 30 天）
bk97
bk97 2017-1-25
编辑： Peng Li 2018-1-18
I have t solve this equation (t^2)*y'' - 2*t*y' + 2*y = (t^3)*log(t) to solve first and secondly to compare the solutions with the theoretical solution y(t)= (7/4)*t + (t^3/2)*log(t) - (3/4)*t^3 ( 1<=t<=2, y'(1)=0 , y(1)=1, step h=0.05 ). I have a possible formula in my book but it doesn't seem to work. Is it possible for someone tο write an indicative formula, please?
  4 个评论
显示 2更早的评论
James Tursa
James Tursa 2017-1-25
编辑：James Tursa 2017-1-25
Why are you using symbolic variables? RK4 is a numerical technique used to solve initial value problems numerically. I.e., you would use it to find the values of your variables at future times. So I don't understand your use of symbols. What are you really trying to do here? E.g., seems like you could use the 2nd order ODE example in this link as a guide for solving your problem:
https://www.mathworks.com/help/matlab/ref/ode45.html?searchHighlight=ode45&s_tid=doc_srchtitle
bk97
bk97 2017-1-25
编辑：bk97 2017-1-25
h=0.05;
t=1:h:2;
y=zeros(1,length(t));
dy=diff(y,t);
y(1)=1;
dy(1)==0;
f=@(t) 2*dy+(2*y)/t+t^2*log(t);
for i=1:(length(t)-1)
k1 = f(t(i),y(i));
k2 = f(t(i)+0.5*h,y(i)+0.5*h*k1);
k3 = f((t(i)+0.5*h),(y(i)+0.5*h*k2));
k4 = f((t(i)+h),(y(i)+k3*h));
y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
i have done this and this error apears "Error using diff Difference order N must be a positive integer scalar." can you help me?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

James Tursa
James Tursa 2017-1-26
编辑：James Tursa 2017-1-26
OK, I will offer a bit more help here (well, actually a lot more help). Your most immediate problem is that you are treating your 2nd order ODE problem as if it is a 1st order ODE problem. So all of your stuff involving y(i) and y(i+1) etc is wrong because that is what you would do for a 1st order ODE (the result at each time step is a scalar). You need to look again at the 2nd order ODE example in this link that I have already given to you:
https://www.mathworks.com/help/matlab/ref/ode45.html?searchHighlight=ode45&s_tid=doc_srchtitle
In there you will see that for the 2nd order ODE problem the "y" solution at each step is in fact a 2-element vector. That was the hint that you needed to pull from this example. Your solution state is a 2-element vector, not a scalar. So all of your equations need to be rewritten with this in mind. The setup is as follows:
y is a 2-element vector
y(1) is your original y
y(2) is y'
So, start with your original equation:
(t^2)*y'' - 2*t*y' + 2*y = (t^3)*log(t)
Solve this for y''
y'' = (2*t*y' - 2*y + (t^3)*log(t)) / t^2
or simplified a bit if you want
y'' = 2*y'/t - 2*y/t^2 + t*log(t)
So notionally if you start with
y(1) = y
y(2) = y'
Then the derivative of this (call it dy) will be
dy(1) = y'
dy(2) = y''
Substituting, you get
dy(1) = y(2)
dy(2) = 2*y(2)/t - 2*y(1)/t^2 + t*log(t)
These expressions for the elements of dy are the equivalent of what your "f" function is in your posted code. So you just need to rewrite your "f" function to take a 2-element y vector and time t, and output a 2-element dy vector. Then use that in your RK4 scheme, keeping in mind that at each point the state is a 2-element vector.
So, changing your scalar y to 2-element y and correcting some of your typos gives you the following framework:
h = 0.05;
t = 1:h:2;
y = zeros(2,length(t)); % <-- changed 1 to 2, each column y(:,i) is the state at time t(i)
y(1,1) = 1; % <-- The initial value of y at time 1
y(2,1) = 0; % <-- The initial value of y' at time 1
f = ____________; % <-- Derivative takes time t and a 2-element y and returns a 2-element result
for i=1:(length(t)-1) % At each step in the loop below, changed y(i) to y(:,i) to accomodate 2-element results
k1 = f( t(i) , y(:,i) );
k2 = f( t(i)+0.5*h, y(:,i)+0.5*h*k1);
k3 = f( t(i)+0.5*h, y(:,i)+0.5*h*k2);
k4 = f( t(i)+ h, y(:,i)+ h*k3);
y(:,i+1) = y(:,i) + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h;
end
Y = (7/4)*t + (t.^3/2).*log(t) - (3/4)*t.^3; % <-- The analytic solution
figure
hold on
plot(t,Y,'b','LineWidth',5); % Plot the analytic solution in thick blue
plot(t,y(1,:),'r','LineWidth',2); % Plot the RK4 solution in thin red
grid on
legend('analytic','RK4');
xlabel('Time')
ylabel('y')
title('Solution to (t^2)*y'''' - 2*t*y'' + 2*y = (t^3)*log(t)');
The only thing you have left to do is fill in the code for the function f. I have given you what you need for this above.
  5 个评论
显示 3更早的评论
James Tursa
James Tursa 2017-1-26
编辑：James Tursa 2017-1-26
dy as a column vector ... remember you are passing in y and returning the derivative of y, which is the dy I spelled out for you above.
Peng Li
Peng Li 2018-1-18
编辑：Peng Li 2018-1-18
Define f as a common function (in another .m file) instead of an anonymous function, which is much easier to implement, as follows. Then the result for the code of @James Tursa can be obtained.
function dy = f( t,y )
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = 2*y(2)./t - 2*y(1)./(t.^2) + t.*log(t);
end
%

请先登录，再进行评论。

Jan
Jan 2017-1-25
A Runge-Kutta as solver for a symbolic function? This is strange. Then your function depends on the inputs y and t, but inside your runge-Kutta-code you call it as f(x) only.
Start with transforming the 2nd order ODE to a set of equations in 1st order. Then omit the "syms", but create the solution numerically. You wil find many working examples when you search for "Matlab runge kutta".
  5 个评论
显示 3更早的评论
bk97
bk97 2017-1-25
编辑：bk97 2017-1-25
h=0.05;
t=1:h:2;
y=zeros(1,length(t));
dy=diff(y,t);
y(1)=1;
dy(1)==0;
f=@(t) 2*dy+(2*y)/t+t^2*log(t);
for i=1:(length(t)-1)
k1 = f(t(i),y(i));
k2 = f(t(i)+0.5*h,y(i)+0.5*h*k1);
k3 = f((t(i)+0.5*h),(y(i)+0.5*h*k2));
k4 = f((t(i)+h),(y(i)+k3*h));
y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
i have done this and this error apears "Error using diff Difference order N must be a positive integer scalar." can you help me?
Jan
Jan 2017-1-26
@bill: The error message means, that in diff the 2nd input means the n'th difference, see doc diff. Therefore diff(y,t) is not meaningful here.
Start with defining a system of 1st order ODEs at first. See e.g. https://www.mathworks.com/help/symbolic/odetovectorfield.html or 2nd order ode to 1st order. Use a function instead of an anonymous function, because it might look easier.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by