《常微分方程的初值问题_Matlab代码》由会员分享，可在线阅读，更多相关《常微分方程的初值问题_Matlab代码（11页珍藏版）》请在金锄头文库上搜索。
1、常微分方程的初值问题常微分方程的初值问题 1.DEEuler 用欧拉法求一阶常微分方程的数值解用欧拉法求一阶常微分方程的数值解 function y = DEEuler(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); x = a:h:b; y(1) = y0; for i=2:N+1 y(i) = y(i-1)+h*Funval(f,varvec,x(i-1), y(i-1); end format short; 2.DEimpEuler 用隐式欧拉法求一阶常微分方程的数值解用隐式欧拉法求一阶常微分方程的数值解 f
2、unction y = DEimpEuler(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 fx = Funval(f,var(1),x(i); gx = y(i-1)+h*fx - varvec(2); y(i) = NewtonRoot(gx,-10,10,eps); end format short; 3.DEModifEuler 用改进欧拉法求一阶常微分方程的数值解用改进欧拉法求一阶常微分方程的数值解
3、function y = DEModifEuler(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 v1 = Funval(f,varvec,x(i-1) y(i-1); t = y(i-1) + h*v1; v2 = Funval(f,varvec,x(i) t); y(i) = y(i-1)+h*(v1+v2)/2; end format short; 4.DELGKT2_mid 用中点法求一阶常微分方程的数
4、值解用中点法求一阶常微分方程的数值解 function y = DELGKT2_mid(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 v1 = Funval(f,varvec,x(i-1) y(i-1); t = y(i-1) + h*v1/2; v2 = Funval(f,varvec,x(i)+h/2 t); y(i) = y(i-1)+h*v2; end format short; 5.DELGKT2_s
5、uen 用休恩法求一阶常微分方程的数值解用休恩法求一阶常微分方程的数值解 function y = DELGKT2_suen(f, h,a,b,y0,varvec) format long; N = uint16(b-a)/h); y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+2*h/3 y(i-1)+K1*2*h/3); y(i) = y(i-1)+h*(K1+3*K2)/
6、4; end format short; 6.DELGKT3_suen 用休恩三阶法求一阶常微分方程的数值解用休恩三阶法求一阶常微分方程的数值解 function y = DELGKT3_suen(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/3 y(i-1)+K1*h/3); K
7、3 = Funval(f,varvec,x(i-1)+2*h/3 y(i-1)+K2*2*h/3); y(i) = y(i-1)+h*(K1+3*K3)/4; end format short; 7.DELGKT3_kuta 用库塔三阶法求一阶常微分方程的数值解用库塔三阶法求一阶常微分方程的数值解 function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Fun
8、val(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2); K3 = Funval(f,varvec,x(i-1)+h y(i-1)-h*K1+K2*2*h); y(i) = y(i-1)+h*(K1+4*K2+K3)/6; end format short; 8.DELGKT4_lungkuta 用经典龙格用经典龙格-库塔法求一阶常微分方程的数值解库塔法求一阶常微分方程的数值解 function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec) format long;
9、 N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2); K3 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K2*h/2); K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*K3); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; end
10、 format short; 9.DELGKT4_jer 用基尔法求一阶常微分方程的数值解用基尔法求一阶常微分方程的数值解 function y = DELGKT4_jer(f, h,a,b,y0,varvec) format long; N = (b-a)/h; C2 = sqrt(2); y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2);
11、 K3 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h*(C2-1)/2+K2*h*(2-C2)/2); K4 = Funval(f,varvec,x(i-1)+h y(i-1)-K2*h*C2/2+K3*h*(2+C2)/2); y(i) = y(i-1)+h*(K1+(2-C2)*K2+(2+C2)*K3+K4)/6; end format short; 10.DELGKT4_qt 用变形龙格用变形龙格-库塔法求一阶常微分方程的数值解库塔法求一阶常微分方程的数值解 function y = DELGKT4_qt(f, h,a,b,y0,varvec) fo
12、rmat long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/3 y(i-1)+K1*h/3); K3 = Funval(f,varvec,x(i-1)+2*h/3 y(i-1)-K1*h/3+K2*h); K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*(K1-K2+K3); y(i) = y(i-1)+h*
13、(K1+3*K2+3*K3+K4)/8; end format short; 11.DELSBRK 用罗赛布诺克半隐式法求一阶常微分方程的数值解用罗赛布诺克半隐式法求一阶常微分方程的数值解 function y = DELSBRK(f, h,a,b,y0,varvec) format long; a1 = 1.40824829; a2 = 0.59175171; c1 = 0.17378667; c2 = c1; w1 = -0.41315432; w2 = 1.41315432; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var
14、 = findsym(f); dy = diff(f, varvec(2); for i=2:N+1 f1 = Funval(f,varvec,x(i-1) y(i-1); dy1 = Funval(dy,varvec,x(i-1) y(i-1); k1 = h*f1/(1-h*a1*dy1); dy2 = Funval(dy,varvec,x(i-1)+c1*h y(i-1)+c2*k1); f2 = Funval(f,varvec,x(i-1)+c1*h y(i-1)+c2*k1); k2 = h*f2/(1-h*a2*dy2); y(i) = y(i-1)+w1*k1+w2*k2; en
15、d format short; 12.DEMS 用默森单步法求一阶常微分方程的数值解用默森单步法求一阶常微分方程的数值解 function y = DEMS(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 fy = Funval(f,varvec,x(i-1) y(i-1); z1 = y(i-1)+h*fy/3; fy1 = Funval(f,varvec,x(i-1)+h/3 z1); z2 = z1+h*(
16、fy1-fy)/6; fy2 = Funval(f,varvec,x(i-1)+h/3 z2); z3 = z2+3*h*(fy2 - 4*fy1/9-fy/9)/8; fy3 = Funval(f,varvec ,x(i-1)+h/2 z3); z4 = z3+2*h*(fy3- 15*fy2/16 + 3*fy/16); fy4 = Funval(f,varvec ,x(i-1)+h z4); y(i) = z4+h*(fy4 - 8*fy3 + 9*fy2 - 2*fy)/6; end format short; 13.DEMiren 用米尔恩法求一阶常微分方程的数值解用米尔恩法求一阶常微分方程的数值解 function y = DEMiren(f, h,a,b,y0,varvec) format long; N = (b-a)/h; y = zeros(N+1,1); x = a:h
