四階龍格-庫(kù)塔法求解常微分方程的初值問(wèn)題-matlab通用程序

1、精選優(yōu)質(zhì)文檔-傾情為你奉上參考教材數(shù)值分析李乃成.梅立泉clearclcformat longm=input('請(qǐng)輸入常微分方程的階數(shù)m=');a=input('請(qǐng)輸入x下限a=');b=input('請(qǐng)輸入x上限b=');h=input('請(qǐng)輸入步長(zhǎng)h=');ym=input('令y(1,1)=y,y（2,1）=y,y(3,1)=y.請(qǐng)輸入ym=','s');     %輸入的時(shí)候必須按照這個(gè)形式輸入y1=y(1,1);if m=1       &

2、#160;                                       %一階初值問(wèn)題單獨(dú)求解    mm=(b-a)/h;    y(1,1)=input('請(qǐng)輸入在初值點(diǎn)的函數(shù)值f(a)=');    x=a;    y11(1)=y(1,1);    for k1=2:(

3、mm+1)        y1=y(1,1);        K(1,1)=h*(eval(ym);                         %計(jì)算K1        x=x+h/2;        y(1,1)=y1+K(1,1)/2;        y1=y(1,1)

4、;        K(1,2)=h*(eval(ym);                         %計(jì)算K2        x=x;        y(1,1)=y1+K(1,2)/2-K(1,1)/2;        y1=y(1,1);        K(1,3)=h

5、*(eval(ym);                          %計(jì)算K3        x=x+h/2;        y(1,1)=y1+K(1,3)-K(1,2)/2;        y1=y(1,1);        K(1,4)=h*(eval(ym);    

6、60;                    %計(jì)算K4        y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4)/6;        y(1,1)=y11(k1);        x=a+(k1-1)*h;           endy11else  

7、;                                               %高階初值問(wèn)題    mm=(b-a)/h;                              

8、       %一共要求解mm個(gè)數(shù)據(jù)點(diǎn)    for k2=1:m                                      %讀取初值條件        fprintf('請(qǐng)輸入%d階導(dǎo)數(shù)的初值f(%d)(a)=n',(k2-1),(k2-1);     

9、   y(k2,1)=input('=');    end    for k2=1:m                                                 y22(1,k2)=y(k2,1);        

10、60;                %先把初值保存在矩陣y22(m,n)中,m表示第幾個(gè)所求點(diǎn)，n表示第n階初值    end    x=a;    for k4=2:(mm+1)                                   %求解mm個(gè)數(shù)據(jù)點(diǎn)的循環(huán) 

11、      for k=1:(m-1)                                 %計(jì)算K1，包括每一階的K1            K(k,1)=h*y(k+1,1);                       &

12、#160;%y(k+1,1)中k+1表示第k+1階，1表示第一個(gè)點(diǎn)；K(k,1)中k表示階數(shù)，1表示K1        end        K(m,1)=h*(eval(ym);        x=x+h/2;                                      %

13、求解K1之前，先重新對(duì)x和y賦值        for k3=1:m                          y(k3,1)=y(k3,1)+K(k3,1)/2;        end        for k=1:(m-1)                

14、                 %計(jì)算K2            K(k,2)=h*y(k+1,1);        end        K(m,2)=h*(eval(ym);        x=x;        for k3=1:m      

15、0;     y(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;        end        for k=1:(m-1)                                  %計(jì)算K3            K(k,3)=h*

16、y(k+1,1);        end        K(m,3)=h*(eval(ym);        x=x+h/2;        for k3=1:m            y(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2;  %這里容易出錯(cuò)        end    

17、0;   for k=1:(m-1)                                     %計(jì)算K4            K(k,4)=h*y(k+1,1);        end        K(m,4)=h*(eval(ym); 

18、      for k5=1:m            y22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4)/6;         %這里，除了要求出下一個(gè)點(diǎn)的數(shù)值，還要求出相應(yīng)的導(dǎo)數(shù)值        end        for k6=1:m                                       %除了對(duì)y(1,1)重新賦值外，還要對(duì)y(2,1)等重新賦值            y(k6,1)=y22(k4,k6);         en