常微分方程組初值問題數(shù)值解的實現(xiàn)和算法分析(編輯修改稿)
2025-02-07 03:32 本頁面
【文章內(nèi)容簡介】 步數(shù) if nargin5 n=100。 end r=size(y0)。 r=r(1)。 s=size(x0)。 s=s(1)。 r=r+s。 T=zeros(r,n+1)。 T(:,1)=[y0。x0]。 for t=2:n+1 k1=feval(f,[T(1:rs,t1)。x0])。 k2=feval(f,[k1*(h/2)+T(1:rs,t1)。x0+h/2])。 k3=feval(f,[k2*(h/2)+T(1:rs,t1)。x0+h/2])。 k4=feval(f,[k3*h+T(1:rs,t1)。x0+h])。 x0=x0+h。 T(:,t)=[T(1:rs,t1)+(k1+k2*2+k3*2+k4)*(h/6)。x0]。 End [2] Euler 改進方法 Matlab 編程 實現(xiàn) 用 Euler改進方法編寫 Matlab程序為 ： function [x,y]=eulerpro(fun,x0,xfinal,y0,n)。 h=(xfinalx0)/n。 x(1)=x0。y(1)=y0。 for i=1:n x(i+1)=x(i)+h。 課程設(shè)計說明書 （論文） 第 6 頁 y1=y(i)+h*feval(fun,x(i),y(i))。 y2=y(i)+h*feval(fun,x(i+1),y1)。 y(i+1)=(y1+y2)/2。 end 5 編程解決 輸入計算題目 function dY=dydx(X,Y) dY(1)=*Y(1)1000*Y(1)*Y(2)。 dY(2)=2500*Y(2)*Y(3)。 dY(3)=*Y(1)1000*Y(1)*Y(2)2500*Y(2)*Y(3) dY=[dY(1)。dY(2)。dY(3)]。 用 RungeKutta 方法的 Matlab 編程解法 function [k,X,Y,wucha,P]=RK4z(dydx,a,b,CT,h) n=fix((ba)/h)。 X=zeros(n+1,1)。 Y=zeros(n+1,length(CT))。 X=a:h:b。 Y(1,:)= CT39。 for k=1:n k1=feval(dydx,X(k),Y(k,:)) x2=X(k)+h/2。y2=Y(k,:)39。+k1*h/2。 k2=feval(dydx,x2,y2)。 k3=feval(dydx,x2,Y(k,:)39。+k2*h/2)。 k4=feval(dydx, X(k)+h,Y(k,:)39。+k3*h)。 Y(k+1,:)=Y(k,:)+h*(k139。+2*k239。+2*k339。+k439。)/6。k=k+1。 end for k=2:n+1 wucha(k)=norm(Y(k)Y(k1))。 k=k+1。 課程設(shè)計說明書 （論文） 第 7 頁 end X=X(1:n+1)。Y=Y(1:n+1,:)。k=1:n+1。wucha=wucha(1:k,:)。 P=[k39。,X39。,Y,wucha39。]。 調(diào)用 ，在 MATLAB工作窗口輸入程序 CT=[1。1。0]。h=。 [k,X,Y,wucha,P]=RK4z(@dydx,0,CT,h), H=[,]。 [x,y]=ode15s(39。dydx39。,H,CT)。 plot(X,Y(:,1),39。g39。,x,y(:,1),39。bo39。,X,Y(:,2),39。m:39。,x,y(:,2),39。cp39。,X,Y(:,3),39。r.39。,x,y(:,3),39。kd39。) xlabel(39。軸 \it x39。)。 ylabel(39。軸 \it y39。) title(39。分別用自定義函數(shù)和 ode15s 函數(shù)求解剛性方程方程組的圖形 39。) legend(39。用 RK4z 函數(shù)解剛性方程的 y1 的曲線 39。,39。用 ode15s 函數(shù)解剛性方程的 y1 的曲線 39。,39。用 RK4z 函數(shù)解剛性方程的 y2 的曲線 39。,39。用 ode15s 函數(shù)解剛性方程的 y2 的曲線 39。, 39。用 RK4z 函數(shù)解剛性方程的 y3 的曲線 39。,39。用 ode15s 函數(shù)解剛性方程的 y3 的曲線 39。) 用改進 Euler 方法的 Matlab 編程解法 function [x,y]=eulerpro(dydx,a,b,CT,h)。 n=fix((ba)/h)。 x=zeros(n+1,1)。 y=zeros(n+1,length(CT))。 y(1,:)= CT39。 for i=1:n x(i+1)=x(i)+h。 f1=feval(dydx,x(i),y(i,:))。 y1(i+1,:)=y(i,:)+h*f139。 f2=feval(dydx,x(i+1),y1(i+1,:))。 y2(i+1,:)=y(i,:)+h*f239。 y(i+1,:)=(y1(i+1,:)+y2(i+1,:))/2。 課程設(shè)計說明書 （論文） 第 8 頁 end 調(diào)用 ，在 MATLAB工作窗口輸入程序 CT=[1。1。0]。h=。a=0。b=。 [X,Y]=eulerpro(@dydx,a,b,CT,h), xlabel(39。軸 \it x39。)。 ylabel(39。軸 \it y39。) H=[,]。 [x,y]=ode15s(39。dy12339。,H,CT)。 plot(X,Y(:,1),39。g39。,x,y(:,1),39。bo39。,X,Y(:,2),39。m:39。,x,y(:,2),39。cp39。,X,Y(:,3),39。r.39。,x,y(:,3),39。kd39。) xlabel(39。軸 \it x39。)
試題試卷相關(guān)推薦
鄂ICP備17016276號-1