學院：數學與資訊科學
專業：資訊與計算科學
姓名：鄭思義
學號：201216524
課程：常微分方程數值解
實驗一：常微分方程的數值解法
1、分別用euler法、改進的euler法（預報校正格式）和s—k法求解初值問題。（h=0.1）並與真解作比較。
1.1實驗**：
%尤拉法
function [x,y]=naeuler(dyfun,xspan,y0,h)
%dyfun是常微分方程，xspan是x的取值範圍，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%改進的尤拉法
function [x,m,y]=naeuler2(dyfun,xspan,y0,h)
%dyfun是常微分方程，xspan是x的取值範圍，y0是初值，h是步長。
%返回值x為x取值，m為預報解，y為校正解
x=xspan(1):h:xspan(2);
y(1)=y0;
m=zeros(length(x)-1,1);
for n=1:length(x)-1
k1=feval(dyfun,x(n),y(n));
y(n+1)=y(n)+h*k1;
m(n)=y(n+1);
k2=feval(dyfun,x(n+1),y(n+1));
y(n+1)=y(n)+h*(k1+k2)/2;
end%四階s—k法
function [x,y]=rk(dyfun,xspan,y0,h)
%dyfun是常微分方程，xspan是x的取值範圍，y0是初值，h是步長。
x=xspan(1):h:xspan(2);
y(1)=y0;
for n=1:length(x)-1
k1=feval(dyfun,x(n),y(n));
k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2);
k3=feval(dyfun,x(n)+h/2,y(n)+(h*k2)/2);
k4=feval(dyfun,x(n)+h,y(n)+h*k3);
y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4
end%主程式
x=[0:0.1:1];y=exp(-x)+x;
dyfun=inline('-y+x+1');
[x1,y1]=naeuler(dyfun,[0,1],1,0.1);
[x2,m,y2]=naeuler2(dyfun,[0,1],1,0.1);
[x3,y3]=rk(dyfun,[0,1],1,0.1);
plot(x,y,'r',x1,y1,'+',x2,y2,'*',x3,y3,'o');
xlabel('x');ylabel('y');
legend('y為真解','y1為尤拉解','y2為改進尤拉解','y3為s—k解','location','northwest');
1.2實驗結果：
2、選取一種理論上收斂但是不穩定的演算法對問題1進行計算，並與真解作比較。（選改進的尤拉法）
2.1實驗思路：演算法的穩定性是與步長h密切相關的。
而對於問題一而言，取定步長h=0.1不論是單步法或低階多步法都是穩定的演算法。所以考慮改變h取值範圍，藉此分析不同步長會對結果造成什麼影響。
故依次採用h=2.0、2.2、2.
4、2.6的改進尤拉法。
2.2實驗**：
%%主程式
x=[0:3:30];y=exp(-x)+x;
dyfun=inline('-y+x+1');
[x1,m1,y1]=naeuler2(dyfun,[0,20],1,2);
[x2,m2,y2]=naeuler2(dyfun,[0,22],1,2.2);
[x3,m3,y3]=naeuler2(dyfun,[0,24],1,2.4);
[x4,m4,y4]=naeuler2(dyfun,[0,26],1,2.6);
subplot(2,2,1)
plot(x,y,'r',x1,y1,'+');xlabel('h=2.0');
subplot(2,2,2)
plot(x,y,'r',x2,y2,'+');xlabel('h=2.2');
subplot(2,2,3)
plot(x,y,'r',x3,y3,'+');xlabel('h=2.4');
subplot(2,2,4)
plot(x,y,'r',x4,y4,'+');xlabel('h=2.6');
2.3實驗結果：
實驗結果分析：
從實驗1結果可以看出，在演算法滿足收斂性和穩定性的前提下，eluer法雖然計算並不複雜，凡是精度不足，反觀改進的eluer法和s—k法雖然計算略微複雜但是結果很精確。
實驗2改變了步長，導致演算法理論上收斂但是不滿足穩定性。結果表示步長h越大，結果越失真。對於同乙個問題，步長h的選取變得尤為重要，這三種單步法的絕對穩定區間並不一樣，所以並沒有一種方法是萬能的，我們應該根據不同的步長來選取合適的方法。
實驗二：ritz-galerkin方法與有限差分法
1、用中心差分格式求解邊值問題
取步長h=0.1，並與真解作比較。
1.1實驗**：
%中心差分法
function u=fdm(xspan,y0,y1,h)
%xspan為x取值範圍，y0,y1為邊界條件，h為步長
n=1/h;
d=zeros(1,n-1);
for i=1:n
x(i)=xspan(1)+i*h;
q(i)=1;
f(i)=x(i);
endfor i=1:n-1
d(i)=q(i)*h*h+2;
end a=diag(d);
b=zeros(n-1);
c=zeros(n-1);
for i=1:n-2
b(i+1,i)=-1;
endfor i=1:n-2
c(i,i+1)=-1;
enda=a+b+c;
for i=2:n-2
b(i,1)=f(i)*h*h;
endb(1,1)=f(1)*h*h+y0;
b(n-1,1)=f(n-1)*h*h+y1;
u= inv(a)*b;
%主程式
x=0:0.1:1;
y=x+(exp(1)*exp(-x))/(exp(2)-1)-(exp(1)*exp(x))/(exp(2)-1);
y1=fdm([0,1],0,0,0.1);
y1=[0,y1',0];
plot(x,y,'r',x,y1,'+')
xlabel('x');ylabel('y');
legend('y真解','y1中心差分法','location','northwest');
1.2實驗結果：
2、用ritz-galerkin方法求解上述問題，並與真值作比較，列表畫圖。
2.1實驗**：
%ritz_galerkin法
function vu=ritz_galerkin(x0,y0,x1,y1,h)
%x0,x1為x取值範圍，y0,y1為邊界條件，h為步長
n=1/h;
syms x;
for i=1:n
fai(i)=x*(1-x)*(x^(i-1));
dfai(i)=diff(x*(1-x)*(x^(i-1)));
endfor i=1:n
for j=1:n
fun=dfai(i)*dfai(j)+fai(i)*fai(j);
a(i,j)=int(fun,x,0,1);
endfun=x*fai(i)+dfai(i);
f(i)=int(fun,x,0,1);
endc=inv(a)*f';
product=c.*fai';
sum=0;
for i=1:n
sum=sum+product(i);
endvu=;
for y=0:h:1
v=subs(sum,x,y);
vu=[vu,v];
endy=0:h:1;
yy=0:0.1:1;
u=sin(yy)/sin(1)-yy;
常微分方程數值解實驗報告
2012年10月17日 實驗目的 1 通過用matlab程式設計運用euler方法及其改進方法求解常微分方程初值問題，更進一步掌握常微分方程及其數值解法課程的理論內容，加深對數值解法的理解。2 熟悉matlab程式設計環境。實驗內容 1.實驗題目 運用euler方法及其改進方法求解常微分方程初值問題...
MATLAB實驗報告 常微分方程數值解
專業序號姓名日期 實驗3 常微分方程數值解 實驗目的 1 掌握用matlab求微分方程初值問題數值解的方法 2 通過例項學習微分方程模型解決簡化的實際問題 3 了解尤拉方法和龍格庫塔方法的基本思想。實驗內容 用尤拉方法和龍格庫塔方法求下列微分方程初值問題的數值解，畫出解的圖形，對結果進行分析比較 解...
微分方程數值解實驗
課程設計報告 班級姓名學號 成績 2017年 6月 21 日 目錄一 摘要 1 二 常微分方程數值解 2 2.1 4階runge kutta法和adams4階外插法的基本思路 2 2.2 演算法流程圖 2 2.3 用matlab編寫源程式 2 2.4 常微分方程數值解法應用舉例 4 三 常係數擴散方...