光电工程师社区

标题: 求教龙格库塔法解速率方程的实例 [打印本页]

作者: 双水 时间: 2009-4-15 13:34
标题: 求教龙格库塔法解速率方程的实例
RT求教龙格库塔法解速率方程的实例,MATLAB实现。

作者: teapot 时间: 2009-4-15 15:09
这个不是有ode的系列命令吗？直接用就行了啊
难道lz要自编的算法程序？

作者: linjipeng 时间: 2009-4-16 10:38
用ｍａｔｌａｂ的帮助ｏｄｅ４５，可以查到实例。误差个人觉得不用很在意，使用默认值就可以

作者: jieweili 时间: 2009-4-16 11:47
%四阶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
??? Undefined function or variable 'N'.

作者: 泪光 时间: 2012-9-19 09:08
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

作者: kcatchy 时间: 2012-11-16 22:08
其实如果不用ode45自己编写也很快，单步循环几十次就可以解决的。

作者: laisai0623 时间: 2012-11-21 19:56
个人认为不用ode45自己编写也好很快

欢迎光临光电工程师社区 (http://bbs.oecr.com/)