求教一个关于RK计算二阶非线性微分方程的问题

*午夜流星*

自己变了一个4阶的RK程序，在求解一个二阶非线性方程（如下）时遇到点问题，不知道是不是迭代次数的问题，就是求解不出来，可是用ode45就可以得到很好的效果。自己的程序以及微分方程文件请见下面，请高手赐教。

function [x,y] = 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)+[0:N]'*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=[y(2);(-1-C*y(2)/y(1)-2*D/y(1)-1.5*y(2)^2)/y(1)];

TBE_Legend

自己变了一个4阶的RK程序，在求解一个二阶非线性方程（如下）时遇到点问题，不知道是不是迭代次数的问题，就是求解不出来，可是用ode45就可以得到很好的效果。自己的程序以及微分方程文件请见下面，请高手赐教。
maple code：

1. restart;
   
2. ode:=b(t)*diff(b(t),t,t)+3/2*diff(b(t),t)^2=-1;
   
3. nsol := dsolve({ode, b(0) = 1,D[1](b)(0)=0}, b(t), numeric);
   
4. with(plots):
   
5. odeplot(nsol,[[t,b(t)],[t,diff(b(t),t)],[t,diff(b(t),t,t)]], 0 .. 1,thickness=2,axes=boxed,legend=["b","b'","b''"]);
   
6. 
   
7. dsolve(b(t)*(diff(diff(b(t), t), t))+(3/2)*(diff(b(t), t))^2 = -1, b(t));
   
8. 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

2阶方程只需要两个个初始条件。

扔掉beta'''(0)=0,用其他两个初始条件接，mmtc和maple结果都说：该方程在t=0.9...出无解，即改点是奇异点。

第二个图是隐式的解析解和6阶级数解

不会用matlab，不知道你ode45怎么求解的，能贴下代码吗？

TBE_Legend

2阶方程只需要两个个初始条件。

扔掉beta'''(0)=0,用其他两个初始条件接，mmtc和maple结果都说：该方程在t=0.9...出无解，即改点是奇异点。

第二个图是隐式的解析解和6阶级数解

不会用matlab，不知道你 ...
TBE_Legend 发表于 2009-12-22 21:50

照猫画虎，用了下matlab。

1. clear all
   
2. clc
   
3. M=[1 0 0;0 1 0;0 0 0]
   
4. y0=[1,0,0]
   
5. tspan=[0,0.9]
   
6. options=odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6],'Vectorized','on')
   
7. [t,y]=ode15s(@f,tspan,y0,options)
   
8. 
   
9. subplot(3,1,1),plot(t,y(:,1),'-r','LineWidth',1.5),ylabel('b'),grid on;
   
10. subplot(3,1,2);plot(t,y(:,2),'-b','LineWidth',1.5);ylabel('db'),grid on;
    
11. subplot(3,1,3);plot(t,y(:,3),'-m','LineWidth',1.5);ylabel('dbb'),grid on;
    
12. 
    
13. % ------------ function f -----------------
    
14. function out=f(t,y);
    
15. out=[y(2,:)
    
16.      y(3,:)
    
17.      y(1,:).*y(3,:)+3/2*y(2,:).^2+1];
    
18. %------------------------------------------
    

复制代码

*午夜流星*

真是多谢斑竹，我再学习学习
3# TBE_Legend

*午夜流星*

ode45的求解程序如下：

>> [x,y]=ode45('bubble_contraction',[0 1],[1 0]);
>> 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

oyanglove5212

多谢楼主的问题，多谢斑竹的回答，把ode15s好好学了下

*午夜流星*

还发现一个很有意思的事情，就是ode45的求解范围永远都比我的RK广，比如C=4;
D=0.025;我的RK只能求解到t[0 9]，而ode45可以求接到t[0 10]，呵呵不知道怎么回事，运行结果如下：
>> [x,y]=PDE_RK_4('bubble_contraction',[0 9],[1 0]);
>> subplot(3,1,1),plot(x,y(:,1),'r'),title('RK,t [0 9]')
>> [x,y]=PDE_RK_4('bubble_contraction',[0 10],[1 0]);
>> subplot(3,1,2),plot(x,y(:,1),'k'),title('RK,t [0 10]')
>> [x,y]=ode45('bubble_contraction',[0 10],[1 0]);
>> subplot(3,1,3),plot(x,y(:,1),'k'),title('ODE45,t [0 10]')




3# TBE_Legend

messenger

看了一下，你编的程序应该没有错误。得不到正确解，应该是RK4算法误差积累的问题。Matlab求解是自适应步长，所以没有这个问题。

另外，将beta除到下面，方程变形为(-1-1.5*y(2)^2)/y(1)本身就不太合理，这使得方程求解多出来一个奇点。

*午夜流星*

对于方程的变形，请问斑竹能给指点一下吗，呵呵
谢谢了
8# messenger