求教龙格库塔法解速率方程的实例

双水 · 发表于 2009-4-15 13:34:42

RT求教龙格库塔法解速率方程的实例,MATLAB实现。

teapot · 发表于 2009-4-15 15:09:09

这个不是有ode的系列命令吗？直接用就行了啊
难道lz要自编的算法程序？

linjipeng · 发表于 2009-4-16 10:38:58

用ｍａｔｌａｂ的帮助ｏｄｅ４５，可以查到实例。误差个人觉得不用很在意，使用默认值就可以

jieweili · 发表于 2009-4-16 11:47:36

%四阶Runge_kutta算法
x=zeros(1,10001);
for i=1:N
K1=h*[a*x(i)-b*x(i)^3+x3(i+1)];
K2=h*[a*(x(i)+K1/2)-b*(x(i)+K1/2)^3+x3(i+1)];
K3=h*[a*(x(i)+K2/2)-b*(x(i)+K2/2)^3+x3(i+1)];
K4=h*[a*(x(i)+K3)-b*(x(i)+K3)^3+x3(i+1)];
x(i+1)=x(i)+K1/6+K2/3+K3/3+K4/6;
end

istanbul · 发表于 2010-5-15 16:34:27

??? Undefined function or variable 'N'.

泪光 · 发表于 2012-9-19 09:08:11

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
调用的子函数以及其调用语句：
function dy=test_fun(x,y)
dy = zeros(3,1);%初始化列向量
dy(1) = y(2) * y(3);
dy(2) = -y(1) + y(3);
dy(3) = -0.51 * y(1) * y(2);
对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：
[T,F] = ode45(@test_fun,[0 15],[1 1 3]);
subplot(121)
plot(T,F)%Matlab自带的ode45函数效果
title('ode45函数效果')
[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数
subplot(122)
plot(T1,F1)%自编的龙格库塔函数效

hkgd · 发表于 2012-11-12 15:35:31

kcatchy · 发表于 2012-11-16 22:08:19

其实如果不用ode45自己编写也很快，单步循环几十次就可以解决的。

laisai0623 · 发表于 2012-11-21 19:56:24

个人认为不用ode45自己编写也好很快


热搜: 防雾膜棒料红外反射膜光电二极管镀膜机反射膜光谱仪光电探头真空镀银滤光片光器件应用光学 zemax激光远心镜头 sio2 增透膜离子源栅网 zemax 论文高反膜半反半透膜

求教龙格库塔法解速率方程的实例

发表回复