Matlab常微分方程數值解法（1）

實驗目的

用Matlab實現歐拉法、后退歐拉法、梯形方法和改進歐拉公式

實驗要求

1. 給出歐拉法、后退歐拉法、梯形方法和改進歐拉公式算法

2. 用Matlab實現歐拉法、后退歐拉法、梯形方法和改進歐拉公式

實驗內容

 實驗步驟

　　（1）歐拉法算法，

　　

　　MATLAB實現，

%數值解常微分方程歐拉算法
%例子：dyfun=inline('y-2*x/y');[x,y]=euler2(dyfun,[0,1],1,0.2);
%輸入：函數dfun(x,y)，求解區間xspan[x0,xN]，初值y0，步長h
%輸出：節點x，數值解y
function [x,y]=euler2(dyfun,xspan,y0,h)
x=xspan(1):h:xspan(2);y(1)=y0;
for n=1:length(x)-1
    y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));
end
x=x';y=y';
end

　　求解【題目】，

 （2）后退歐拉法算法，

　　MATLAB實現

%數值解常微分后退歐拉法算法
%例子：dfun=inline('x+y','x','y');[x,y]=eulerh1(dfun,0,1,0.02,5)
%輸入：函數dfun(x,y)，初值x0,y0,步長h,維度N
%輸出：結點x和數值解y
function [x,y]=eulerh1(dfun,x0,y0,h,N)
x=zeros(1,N+1);
y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
    x(n+1)=x(n)+h;
    z0=y(n)+h*dfun(x(n),y(n));
    for k=1:3
        z1=y(n)+h*dfun(x(n+1),z0);
        if abs(z1-z0)<1e-6
            break;
        end
        z0=z1;
    end
    y(n+1)=z1;
end

　　求解【題目】，

　　

 （3）梯形方法算法，

 　　MATLAB實現，

%數值解常微分梯形歐拉法算法
%例子：dfun=inline('x+y','x','y');[x,y]=eulert(dfun,0,1,0.02,5)
%輸入：函數dfun(x,y)，初值x0,y0,步長h,維度N
%輸出：結點x和數值解y
function [x,y]=eulert(dfun,x0,y0,h,N)
x=zeros(1,N+1);
y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
    x(n+1)=x(n)+h;
    z0=y(n)+h*dfun(x(n),y(n));
    for k=1:3
        z1=y(n)+(h/2)*(dfun(x(n),y(n))+dfun(x(n+1),z0));
        if abs(z1-z0)<1e-6
            break;
        end
        z0=z1;
    end
    y(n+1)=z1;
end

 　　求解【題目】，

 　　

（2）改進歐拉公式算法，

 　　MATLAB實現，

%數值解常微分改進歐拉法算法
%例子：dfun=inline('x+y','x','y');[x,y]=eulerh1(dfun,0,1,0.02,5)
%輸入：函數dfun(x,y)，初值x0,y0,步長h,維度N
%輸出：結點x和數值解y
function [x,y]=eulerg2(dfun,x0,y0,h,N)
x=zeros(1,N+1);
y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
    x(n+1)=x(n)+h;
    ybar=y(n)+h*dfun(x(n),y(n));
    y(n+1)=y(n)+(h/2)*(dfun(x(n),y(n))+dfun(x(n+1),ybar));
end
end

　　求解【題目】，

　　

求解結果，

　　

 小結

　　就給定的題目並沒有體現出這些算法之間的差異。