如何用改进欧拉方法求解

120397406

求解：
dy/dx=y+2*x/y^2  y(0)=1 x∈[0，2]
怎么用改进欧拉的方法编写程序啊？？？步长要求等于0。1
还要能画出图。。
哪位高手能帮帮忙啊，我要做报告。。

xjzuo

建议看看数值分析的书, 这是很基础的作业.
你也可以搜索一下论坛.

120397406

就是找不到。。。

happy

1. function w = ieulerf(f,t,u,h)
   
2. %
   
3. %  v = ieulerf(f,t,u,h)
   
4. %   
   
5. %  Computes one step of Improved Euler method for solving a first order ODE
   
6. %
   
7. % Input arguments:
   
8. %
   
9. %   f  - string name of function on right hand side of ODE: u' = f(t,u)
   
10. %   t  - current time
    
11. %   u  - current value of solution
    
12. %   h  - step size
    
13. %
    
14. % Output:
    
15. %
    
16. %   v  - updated value of solution
    
17. %
    
18. % See also NUMODE, NEULER, MIDP, IEULER, RK4, NEULERF, MIDPF, RK4F
    
19. %
    
20. 
    
21. g = feval(f,t,u);
    
22. v = u + h*g;
    
23. w = u + h*(g + feval(f,t+h,v))/2;

复制代码

happy

1. function [tz,uz] = ieuler(f,t0,u0,tf,h,c)
   
2. %
   
3. %  [tz,uz] = ieuler(f,t0,u0,tf,h,c)
   
4. %   
   
5. %  Uses Improved Euler method to solve a first order ODE
   
6. %    Plots components of solution in Figure #1
   
7. %    In 2-D case, also plots phase plane trajectory in Figure #2
   
8. %
   
9. % Input arguments:
   
10. %
    
11. %   f  - string name of function on right hand side of ODE: u' = f(t,u)
    
12. %            note: f must return a column vector
    
13. %   t0 - initial time
    
14. %   u0 - initial data u(t0) --- column vector
    
15. %   tf - final time or, if vector, list of times to evaluate solution.
    
16. %   h  - step size
    
17. %
    
18. % Optional arguments:
    
19. %
    
20. %   c  - color for plot
    
21. %
    
22. % Output:
    
23. %
    
24. %   tz - times for evaluation of function as specified by tf, or,
    
25. %        if tf is a scalar, all mesh points
    
26. %   uz - value of solution at times tz.
    
27. %
    
28. % See also NUMODE, NEULER, MIDP, RK4, NEULERF, IEULERF, MIDPF, RK4F
    
29. %
    
30. 
    
31. 
    
32. [tz,uz] = numode('ieulerf',f,t0,u0,tf,h,c);

复制代码

120397406

xnp=0:.1:2;
y21=1;
ynp=[];
for i=1:length(xnp)-1
    y04=y21+diff(xnp([1:2]))*(y21-2*xnp(i)/y21^2);
    ynp=[ynp;y21+.5*diff(xnp([1:2]))*(y21-2*xnp(i)/y21^2+y04-2*xnp(i+1)/y04^2)];
    y21=ynp(end);
end
ynp
我自己写好了。。但是问题是怎么把图象显示出来？谁会啊？？？我输入一直显示        
指数超过矩阵尺寸.