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

2、#160;                                       %一階初值問題單獨求解    mm=(b-a)/h;    y(1,1)=input('請輸入在初值點的函數(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);                         %計算K1        x=x+h/2;        y(1,1)=y1+K(1,1)/2;        y1=y(1,1)

4、;        K(1,2)=h*(eval(ym);                         %計算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);                          %計算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;                    %計算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、;                                               %高階初值問題    mm=(b-a)/h;                              

8、       %一共要求解mm個數(shù)據(jù)點    for k2=1:m                                      %讀取初值條件        fprintf('請輸入%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表示第幾個所求點，n表示第n階初值    end    x=a;    for k4=2:(mm+1)                                   %求解mm個數(shù)據(jù)點的循環(huán) 

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

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

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

14、                 %計算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)                                  %計算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;  %這里容易出錯        end    

17、0;   for k=1:(m-1)                                     %計算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;         %這里，除了要求出下一個點的數(shù)值，還要求出相應(yīng)的導(dǎo)數(shù)值        end        for k6=1:m                                       %除了對y(1,1)重新賦值外，還要對y(2,1)等重新賦值            y(k6,1)=y22(k4,k6);         en