數值分析實驗之常微分方程的數值解法（MATLAB實現）

一、實驗目的

科學技術中常常要求解常微分方程的定解問題，所謂數值解法就是求未知函數在一系列離散點處的近似值。

二、實驗原理

三、實驗程序

  1. 尤拉公式程序

     

 

      2、3、4的尤拉公式的程序參上改寫。

 四、實驗內容

   

 

五、實驗代碼及運行結果

    • MATLAB代碼：

定義函數：
function [A1,A2,B1,B2,C1,C2]=euler23(a,b,n,y0)
%歐拉法解一階常微分方程
%初始條件y0
h = (b-a)/n; %步長h
%區域的左邊界a
%區域的右邊界b
x = a:h:b;
m=length(x);
%前向歐拉法
y = y0;
for i=2:m
    y(i)=y(i-1)+h*oula(x(i-1),y(i-1));
    A1(i)=x(i);
    A2(i)=y(i);
end
plot(x,y,'r-');
hold on;
%改進歐拉法
y = y0;
for i=2:m
    y(i)= y(i-1)+h/2*( oula(x(i-1),y(i-1))+oula(x(i),y(i-1))+h*(oula(x(i-1),x(i-1))));
    B1(i)=x(i);
    B2(i)=y(i);
end
plot(x,y,'m-');
hold on;
%歐拉兩步公式
y=y0;
y(2)=y(1)+h*oula(x(1),y(1));
for i=2:m-1
    y(i+1)=y(i-1)+2*h*oula(x(i),y(i));
    C1(i)=x(i);
    C2(i)=y(i);
end
plot(x,y,'b-');
hold on;
%精確解用作圖
xx = x;
f = dsolve('Dy=-y+1','y(0)=0','x');%求出解析解
y = subs(f,xx); %將xx代入解析解，得到解析解對應的數值
plot(xx,y,'k--');
legend('前向歐拉法','改進歐拉法','歐拉兩步法','解析解');

function f=oula(x,y)
f=-y+1;

命令行窗口：
[A1,A2,B1,B2,C1,C2]=euler23(0,1,10,0)

　　運行結果：

      

  N=50時：

      

 

N=100時：

      

 

   故得精度越大時，幾種方法求解值與准確值越來越接近。

 

     • 另解

clear;
format long;
a = 0;
b = 1;
h = 0.1;
d = 0;
res = forward(a, b, h, d);
x = res(1,:);
y = res(2,:);
xx = x;
f = dsolve('Dy=-y+1','y(0)=0','x');
z  = subs(f,xx); 
y(2,:) = z;
plot(x, y);



function result = forward(a, b, h, y)
    n = (b-a)/h;
    x0 = a;
    x1 = a;
    y0 = y;
    result(1,1) = x0;
    result(2,1) = y0;

    for m = 0:n-1
        x1 = x1 + h;
        f0 = 1-y0;
        d = y0 + h*f0;
        y1 = calculate(y0, x1, d, h);
        %result = calculate(x1, d, h);
        x0 = x1;
        y0 = y1;
        result(1, m+2) = x0;
        result(2, m+2) = y0;

    end
end

function result = calculate(y0, x1, y1, h)

    acc = -6;
    now = 0.0;
    z1 = y1;

    while now >= -6

        z0 = z1;
        f0 =1-z0;
        z1 = y0 + h*f0;
        now = log10(abs(z1-z0));

    end
    result = z1;
end

　　運行結果：