求编程达人帮忙用matlab编程用龙格库塔方法解微分方程

Posted 2023-03-17

用变步长的四阶龙格库塔方法求初值问题 y'= y-2x/y,y(0)=1.的数值解。精度为10的负8次方。并与4阶经典龙格库塔方法（ode45)所得结果比较，并作图。

求程序达人帮忙，另有加分。

给一个程序，朋友试试看吧，呵呵！

功能：用四阶Runge-Kutta 法求解常微分方程
---------------------------------------------
function R=Rungkuta4(f, a, b, n, ya)
% f:微分方程右端函数句柄
% a,b:自变量取值区间的两个端点
% n:区间等分的个数
% ya:函数初值y(a)
% R=[x',y']:自变量X 和解Y 所组成的矩阵
h=(b-a)/n;
x=zeros(,n+1);
y=zeros(1,n+1);
x=a:h:b;
y(1)=ya;
for i=1:n
k1=h*feval(f,x(i),y(i));
k2=h*feval(f,x(i)+h/2,y(i)+k1/2);
k3=h*feval(f,x(i)+h/2,y(i)+k2/2);
k4=h*feval(f,x(i)+h,y(i)+k3);
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6;
end
R=[x',y']; 参考技术A clear
clc
n=4;
h0=(1-0)/n;h=h0;m=0;
y(1)=1;y(n+1)=2;y0(3)=33;n=n/2;
while abs(y(n+1)-y0(3))>((1e-8)/15)
n=n*2;x=0:h:1;
for k=1:n
k1=tfb(x(k),y(k));
k2=tfb((x(k)+h/2),(y(k)+h*k1/2));
k3=tfb((x(k)+h/2),(y(k)+h*k2/2));
k4=tfb((x(k)+h),(y(k)+h*k3));
y(k+1)=y(k)+(h/6)*(k1+2*k2+2*k3+k4);
end
h=h/2;x0(1)=x(n);x0(2)=x(n)+h;y0(1)=y(n);
for k=1:2
k1=tfb(x0(k),y0(k));
k2=tfb((x0(k)+h/2),(y0(k)+h*k1/2));
k3=tfb((x0(k)+h/2),(y0(k)+h*k2/2));
k4=tfb((x0(k)+h),(y0(k)+h*k3));
y0(k+1)=y0(k)+(h/6)*(k1+2*k2+2*k3+k4);
end
m=m+1;
end
[x1 y1]=ode45('tfb',[0 1],[1]);
plot(x,y,'-',x1,y1,'+')
x
y
m
h
x1
y1本回答被提问者采纳