MATLAB: What happens to this programe of Runge-Kutta Method

ode

I want to solve the Initial Value Problem below:
y''+y=0.001x^2,y(0)=0,y'(0)=1.5
Therefor,I let y1=y,and y2=y',it yields a systems of ODE :
y1'=y2
y2'=0.001x^2-y1
Thus,I want to solve this Initial Value Problem by Runge-Kutta Method of 4-th orders with the programe below:
clear;close all;clc
%Initial Condition 
y_0=0;
y_prime_0=1.5;
%empty k vector
k1=[0,0];
k2=[0,0];
k3=[0,0];
k4=[0,0];
%empty y vector
x=zeros(1,102);
y1=zeros(1,102);
y2=zeros(1,102);
y1(1)=y_0;
y2(1)=y_prime_0;
%each step
h=0.5;
%Iteration by Runge-Kutta Method
for i=1:101
    k1(1)=h*y2(i);
    k1(2)=h*(0.001*x(i)^2-y1(i));
    k2(1)=h*(y2(i)+0.5*k1(1));
    k2(2)=h*(0.001*(x(i)+0.5*h)^2-(y1(i)+0.5*k1(2)));
    k3(1)=h*(y2(i)+0.5*k2(1));
    k3(2)=h*(0.001*(x(i)+0.5*h)^2-(y1(i)+0.5*k2(2)));
    k4(1)=h*(y2(i)+k3(1));
    k4(2)=h*(0.001*(x(i)+h)^2-(y1(i)+k3(2)));
    y1(i+1)=y1(i)+(1/6)*(k1(1)+2*k2(1)+2*k3(1)+k4(1));
    y2(i+1)=y2(i)+(1/6)*(k1(2)+2*k2(2)+2*k3(2)+k4(2));
    x(i+1)=x(i)+h;
end
%Numerical Solution
plot(x,y1,'-k');
hold on;
%General Solution
y_c=0.002*cos(x)+1.5*sin(x)+0.001*x.^2-0.002;
plot(x,y_c,'-b');
axis([0,50,-5,50]);
At last,I check the result with general solution y=0.002cosx+1.5sinx+0.001x^2-0.002,but the figure becomes weird:
So what happens to this programe,and how can I fix it?

Best Answer

Welcome to numerical integration! Nothing is wrong with your code. Your step size is too large for your system. Make it smaller, and your results will more closely resemble the general solution. Here I use .05 instead of .5. You will have to increase your number of loops accordingly to compensate (10x here).
%Initial Condition 
y_0=0;
y_prime_0=1.5;
%empty y vector
x=0;
y1(1)=y_0;
y2(1)=y_prime_0;
%each step
h=0.05;
%Iteration by Runge-Kutta Method
for i=1:1001
    k1(1)=h*y2(i);
    k1(2)=h*(0.001*x(i)^2-y1(i));
    k2(1)=h*(y2(i)+0.5*k1(1));
    k2(2)=h*(0.001*(x(i)+0.5*h)^2-(y1(i)+0.5*k1(2)));
    k3(1)=h*(y2(i)+0.5*k2(1));
    k3(2)=h*(0.001*(x(i)+0.5*h)^2-(y1(i)+0.5*k2(2)));
    k4(1)=h*(y2(i)+k3(1));
    k4(2)=h*(0.001*(x(i)+h)^2-(y1(i)+k3(2)));
    y1(i+1)=y1(i)+(1/6)*(k1(1)+2*k2(1)+2*k3(1)+k4(1));
    y2(i+1)=y2(i)+(1/6)*(k1(2)+2*k2(2)+2*k3(2)+k4(2));
    x(i+1)=x(i)+h;
end
%Numerical Solution
plot(x,y1,'-k');
hold on;
%General Solution
y_c=0.002*cos(x)+1.5*sin(x)+0.001*x.^2-0.002;
plot(x,y_c,'-b');
axis([0,50,-5,50]);
●
Related Question