...源代码用四阶龙格库塔法解如下微分方程 y'=y-2x/y(0<x<1),y(0...
发布网友
发布时间：2024-09-29 11:12
共1个回答
热心网友
时间：2024-09-29 12:08
function [x,y]=myrk4(ufunc,y0,h,a,b)
%参数: 函数名称，初始值向量，步长，时间起点，时间终点
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值
%按龙格库塔方法进行求解
for ii=1:n
x(ii+1)=x(ii)+h;
k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);
k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);
k4=ufunc(x(ii)+h,y(:,ii)+h*k3);
y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
end
以下是主程序
% y'=y-2x/y (0<x<1),y(0)=1,步长为h=0.2fun = inline('y-2*x/y');
[t1,f1]=myrk4(fun,1,0.2,0,1);%测试时改变test_fun的函数维数，别忘记改变初始值的维数
subplot(211); plot(t1,f1) %自编函数
title('自编函数求解结果')
%用系统自带函数ode45进行比较
[t,f] = ode45(fun,[0 1],1);
subplot(212); plot(t,f);title('ode45求解结果')
