请各位帮我看个关于精度溢出的问题

wwg.no.19

1. function num_euler_formula2
   
2. clc
   
3. %m=input('m=');
   
4. %c=input('c=');
   
5. %k=input('k=');
   
6. %x0=input('x0=');
   
7. %v0=input('v0=');
   
8. 
   
9. %f0=input('f0=');
   
10. %w=input('w=');
    
11. %delt=input('delt=');
    
12. %wn=sqrt(k/m);
    
13. %z=c/2/m/wn;
    
14. %lan=w/wn;
    
15. %x=x0;
    
16. %xd=v0;
    
17. %T=input('T=');
    
18. m=18.2
    
19. c=1.49
    
20. k=43.8
    
21. x0=1
    
22. v0=1
    
23. f0=44.5
    
24. w=15
    
25. delt=1
    
26. T=70;
    
27. x=x0;
    
28. xd=v0;
    
29. wn=sqrt(k/m);
    
30. z=c/2/m/wn;
    
31. lan=w/wn;
    
32. 
    
33. 
    
34. 
    
35. for t=1:delt:T   
    
36. 
    
37.    
    
38.     xdd=(f0*sin(w*t)-k*x0-c*xd)/m;
    
39.     x0=x0+xd*delt;
    
40.     xd=xd+xdd*delt;
    
41.     x(t)=x0;
    
42.     xxd(t)=xd;
    
43. end
    
44. n=T/delt;
    
45. 
    
46. t=1:delt:T
    
47. x
    
48. 
    
49. 
    
50. 
    
51. subplot(1,2,1);
    
52. plot(t,x);
    
53. grid;
    
54. xlabel('time');
    
55. ylabel('displacement');
    
56. title('time and displacement');
    
57. 
    
58. subplot(1,2,2);
    
59. plot(t,xxd);
    
60. grid;
    
61. xlabel('time');
    
62. ylabel('velocity');
    
63. title('time and velocity');

复制代码

用欧拉法求单自由度谐迫响应，但位移X初始值和最终值差太多，导致用e表示时系统把初始值都认为0了，最终画图时初始阶段画出来的都是一条直线。而且位移应该是衰减的，怎么会递增了