want to solve (d^2 y)/(dx^2 )+dy/dx-6y=0 using 4th order Runge-Kutta method with y(0) = 3 and y’(0) = 1

19 次查看(过去 30 天)
My code is :
clc
clear all
h = 0.5;
x = 1:h:5;
y = zeros(2,length(x)); % y vector declaration
y(1,1) = 3; % y value
y(2,1) = 1; % y' value
f = @(x) (dy^2)/(dx^2 )+dy/dx-6*y;
for i=1:(length(x)-1) % calculating y values
z1 = f( x(i) , y(:,i));
z2 = f( x(i)+0.5*h, y(:,i)+0.5*h*z1);
z3 = f( x(i)+0.5*h, y(:,i)+0.5*h*z2);
z4 = f( x(i)+h, y(:,i)+h*z3);
y(:,i+1) = y(:,i) + (1/6)*(z1 + 2*z2 + 2*z3 + z4)*h;
end
%plotting results
Y = x.^2 - x - 6;
figure
hold on
plot(x,Y,'b'); % analytic solution plot--- blue
plot(x,y(1,:),'r'); % runga kutta solution --- Red
grid on
legend('Analytical solution','RK4');
xlabel('time (T)')
ylabel('y')
title('Runga gutta 4th order');
But My output is :
Too many input arguments.
Error in Untitled3 (line 11)
z1 = f( x(i) , y(:,i));
How can i fix this problem sir?

采纳的回答

James Tursa
James Tursa 2021-2-7
编辑:James Tursa 2021-2-7
Change your function handle from this
f = @(x) (dy^2)/(dx^2 )+dy/dx-6*y;
to this
f = @(x,y) [y(2);-y(2)+6*y(1)];
The last part of that function handle is simply solving this
(d^2 y)/(dx^2 )+dy/dx-6y=0
for the 2nd derivative and then plugging in y(2) and y(1) for the other parts.

更多回答(1 个)

Riccardo Bonomelli
Hi, the problem seems to be in the definition of the function f, you defined f as a single variable function or f(x) while, according to your code, f should be a function of two variables, like f(x,y). Furthermore, the dy and dx inside the function f in your definition must be defined in some way, say dx= something and dy= something. I guess you are trying to express a derivative, but you must define it in some way.
  3 个评论
Riccardo Bonomelli
The definition of the derivatives depends on the equation you are trying to solve, in your case you are solving a second order differential equation which requires a system of two first order differential equations to be solved using 4th order Runge Kutta. I'm no expert on the subject, so I am not comfortable explaining the details of the implementation.

请先登录,再进行评论。

Community Treasure Hunt

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

Start Hunting!

Translated by