function [t,xt]=oderk(f,tI,x0,n)
% 龙格库塔 求常微分方程初值问题数值解法, 本次实验我写的是四阶龙格库塔
% f 导数值 tI x的定义域 x0 初值
if nargin<4
n=1;
end
% 定义域区间为2时,将他分为100个区间
if length(tI)==2
t=linspace(tI(1),tI(2));
else
t=tI;
end
switch n
case 1
% 四阶龙格库塔公式
k1=@(t,x,h)f(t,x);
k2=@(t,x,h)f(t+h/2,x+h/2*k1(t,x,h));
k3=@(t,x,h)f(t+h/2,x+h/2*k2(t,x,h));
k4=@(t,x,h)f(t+h,x+h*k3(t,x,h));
w=@(t,x,h)(k1(t,x,h)+2*k2(t,x,h)+2*k3(t,x,h)+k4(t,x,h))/6;
end
%非等距分割,求导
dt=diff(t);
%定义数组,保存xt
xt=zeros(length(x0),length(t));
%初值,保存零点值
xt(:,1)=x0;
%从t-t-1 开始迭代
for k=1:length(t)-1
xt(:,k+1)=xt(:,k)+dt(k)*w(t(k),xt(:,k),dt(k));
end
end
oderk
这篇博客介绍了如何使用四阶龙格库塔方法求解常微分方程初值问题。代码实现了一个函数`oderk`,该函数根据给定的定义域、初始值和阶数,进行非等距分割,计算并返回数值解。函数灵活适应不同长度的定义域,并在区间为2时默认分为100个小区间。整个过程通过矩阵运算和迭代完成。
摘要由CSDN通过智能技术生成