Question

求助一道matlab的题 谢谢大家了

ralston's 方法,
a2=2/3
a1=1/3, p1=q11=3/4
y i+1=y i + [(1/3)k1+(2/3)k2]h
k1=f(xi,yi)
k2=f(xi+(3/4)h, yi+(3/4)k1h)

Answer 1

function [x,y]=Runge_kutta2(f,a,b,x0)

%2阶Runge_kutta解微分方程

%调用格式同ode45，

%f为微分方程函数，a b为积分区间，x0初值

%h为步长，默认为0.001

h=0.001;

xk=a:h:b;

n=(b-a)/h+1;

a=ones(length(x0),1);

y1(:,1)=x0';

for i=1:n

 x1=xk(i)*a;

k1=f(x1,y1(:,i));

k2=f(x1+(3/4)*h*a,y1(:,i)+(3/4)*k1*h);

y1(:,i+1)=y1(:,i)+((1/3)*k1+(2/3)*k2)*h;

y(i,:)=y1(:,i)';

end

x=xk;

--------------------------------------------------------------

odefun=@(t,x)[-10*x(1)*x(3)+x(2)

    10*x(1)*x(3)-x(2)

    -10*x(1)*x(3)+x(2)-2*x(3)];

[t,y]=ode45(odefun,[0 10],[50 0 40]);

[t1,y1]=Runge_kutta2(odefun,0,10,[50 0 40])

subplot(2,1,1),plot(t,y);legend('a-t','b-t','c-t');title('ode45')

subplot(2,1,2),plot(t1,y1);legend('a-t','b-t','c-t');title('Runge_kutta2')

Answer 2

课件的p8.11的公式得给出来。
因为龙格库塔法的系数可以不唯一。