求教一个关于RK计算二阶非线性微分方程的问题
本帖最后由 *午夜流星* 于 2009-12-22 20:36 编辑自己变了一个4阶的RK程序，在求解一个二阶非线性方程（如下）时遇到点问题，不知道是不是迭代次数的问题，就是求解不出来，可是用ode45就可以得到很好的效果。自己的程序以及微分方程文件请见下面，请高手赐教。
function = PDE_RK_4(f,x_span,y0,N)
if nargin<4 | N<=0
N = 10000;
end
if nargin<3
y0=0;
end
y(1,:)=y0(:)';
h = (x_span(2) - x_span(1))/N;
x = x_span(1)+'*h;
for k=1:N
f1 = h*feval(f,x(k),y(k,:));
f1 = f1(:)';
f2 = h*feval(f,x(k) + h/2,y(k,:) + f1/2);
f2 = f2(:)';
f3 = h*feval(f,x(k) + h/2,y(k,:) + f2/2);
f3 = f3(:)';
f4 = h*feval(f,x(k) + h,y(k,:) + f3);
f4 = f4(:)';
y(k + 1,:) = y(k,:) + (f1 + 2*(f2 + f3) + f4)/6;
end
微分方程程序：
function dy=bubble_contraction(t,y)
C=0;
D=0;
dy=; 本帖最后由 TBE_Legend 于 2009-12-23 09:13 编辑
自己变了一个4阶的RK程序，在求解一个二阶非线性方程（如下）时遇到点问题，不知道是不是迭代次数的问题，就是求解不出来，可是用ode45就可以得到很好的效果。自己的程序以及微分方程文件请见下面，请高手赐教。
maple code：
restart;
ode:=b(t)*diff(b(t),t,t)+3/2*diff(b(t),t)^2=-1;
nsol := dsolve({ode, b(0) = 1,D(b)(0)=0}, b(t), numeric);
with(plots):
odeplot(nsol,[,,], 0 .. 1,thickness=2,axes=boxed,legend=["b","b'","b''"]);
dsolve(b(t)*(diff(diff(b(t), t), t))+(3/2)*(diff(b(t), t))^2 = -1, b(t));
dsolve(b(t)*(diff(diff(b(t), t), t))+(3/2)*(diff(b(t), t))^2 = -1, b(t), type = series);
22 ...
*午夜流星* 发表于 2009-12-22 20:35 http://forum.simwe.com/images/common/back.gif
2阶方程只需要两个个初始条件。
扔掉beta'''(0)=0,用其他两个初始条件接，mmtc和maple结果都说：该方程在t=0.9...出无解，即改点是奇异点。
第二个图是隐式的解析解和6阶级数解
不会用matlab，不知道你ode45怎么求解的，能贴下代码吗？
2阶方程只需要两个个初始条件。
扔掉beta'''(0)=0,用其他两个初始条件接，mmtc和maple结果都说：该方程在t=0.9...出无解，即改点是奇异点。
第二个图是隐式的解析解和6阶级数解
不会用matlab，不知道你 ...
TBE_Legend 发表于 2009-12-22 21:50 http://forum.simwe.com/images/common/back.gif
照猫画虎，用了下matlab。
clear all
clc
M=
y0=
tspan=
options=odeset('Mass',M,'RelTol',1e-4,'AbsTol',,'Vectorized','on')
=ode15s(@f,tspan,y0,options)
subplot(3,1,1),plot(t,y(:,1),'-r','LineWidth',1.5),ylabel('b'),grid on;
subplot(3,1,2);plot(t,y(:,2),'-b','LineWidth',1.5);ylabel('db'),grid on;
subplot(3,1,3);plot(t,y(:,3),'-m','LineWidth',1.5);ylabel('dbb'),grid on;
% ------------ function f -----------------
function out=f(t,y);
out=[y(2,:)
y(3,:)
y(1,:).*y(3,:)+3/2*y(2,:).^2+1];
%------------------------------------------
真是多谢斑竹，我再学习学习
3# TBE_Legend ode45的求解程序如下：
>> =ode45('bubble_contraction',,);
>> subplot(2,1,1),plot(x,y(:,1),'-k','LineWidth',1.5),ylabel('b'),grid on;
>> subplot(2,1,2),plot(x,y(:,2),'-r','LineWidth',1.5),ylabel('d(b)'),grid on;
但是还有一个问题就是，如果原来函数中的C或者D不为零，方程就很难求解了，比如C=1，D=0时，还要再跟版主学习，呵呵
2# TBE_Legend 多谢楼主的问题，多谢斑竹的回答，把ode15s好好学了下 还发现一个很有意思的事情，就是ode45的求解范围永远都比我的RK广，比如C=4;
D=0.025;我的RK只能求解到t，而ode45可以求接到t，呵呵不知道怎么回事，运行结果如下：
>> =PDE_RK_4('bubble_contraction',,);
>> subplot(3,1,1),plot(x,y(:,1),'r'),title('RK,t ')
>> =PDE_RK_4('bubble_contraction',,);
>> subplot(3,1,2),plot(x,y(:,1),'k'),title('RK,t ')
>> =ode45('bubble_contraction',,);
>> subplot(3,1,3),plot(x,y(:,1),'k'),title('ODE45,t ')
3# TBE_Legend 看了一下，你编的程序应该没有错误。得不到正确解，应该是RK4算法误差积累的问题。Matlab求解是自适应步长，所以没有这个问题。
另外，将beta除到下面，方程变形为(-1-1.5*y(2)^2)/y(1)本身就不太合理，这使得方程求解多出来一个奇点。 对于方程的变形，请问斑竹能给指点一下吗，呵呵
谢谢了
8# messenger
页:
[1]