【實(shí)驗(yàn)原理】
一迫吐、改進(jìn)的歐拉法(預(yù)報(bào)校正)
1)預(yù)報(bào)校正法基本思路
改進(jìn)歐拉法先用歐拉法求出預(yù)報(bào)值,再利用梯形公式求出校正值溯祸,局部截?cái)嗾`差比歐拉法低了一階叹俏,可以較大程度地提高計(jì)算精度。
2)預(yù)報(bào)校正法算法描述
二嚷往、四階規(guī)范龍格-庫塔公式
【實(shí)驗(yàn)內(nèi)容】
一葛账、回答下面的問題
1.什么是常微分方程的解析解和數(shù)值解?
2.歐拉左矩形公式和右矩形公式是什么,有什么區(qū)別皮仁。
3.解函數(shù)y(x)的光滑性表現(xiàn)在什么地方籍琳,解函數(shù)y(x)“變化劇烈”如何體現(xiàn),請舉例說明贷祈。
4.寫出一個二階龍格-庫塔公式趋急、四階標(biāo)準(zhǔn)龍格-庫塔公式
二、編程計(jì)算課后習(xí)題中規(guī)定題目势誊,交回實(shí)驗(yàn)報(bào)告與計(jì)算
實(shí)驗(yàn)代碼
%定義函數(shù)
function z=f1(x,y)
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
z=1/(x^2)-y/x;
end
%改進(jìn)歐拉法
clear;clc;
h=0.05;
x=1:h:2; %定義函數(shù)范圍
y(1)=1; %初值
for i=1:20
yp=y(i)+h*f1(x(i),y(i));%改進(jìn)(迭代)公式
yc=y(i)+h*f1(x(i+1),yp);
y(i+1)=(yp+yc)/2;%校正值
end
%-----------輸出--------------------
disp('i=0...10的值為:');
for i=1:2:21
fprintf('%f\n',y(i))
end
%-----------輸出--------------------
%四階RK
clear all
a=1;b=2;h=0.1;%邊界及步長
x0=a;y0=1;%初值
n=(b-a)/h;
N=[0:1:n]';
x=zeros(n+1,1);
y=zeros(n,1);%用矩陣存儲
x(1)=x0;y(1)=y0;
for i=1:n
x(i+1)=x(i)+h;
end
%----------------四階RK--------------------
for i=1:n
k1=h*f1(x(i),y(i));
k2=h*f1(x(i)+h/2,y(i)+k1/2);
k3=h*f1(x(i)+h/2,y(i)+k2/2);
k4=h*f1(x(i)+h,y(i)+k3);
y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6;
end
%----------------四階RK--------------------
[N,x,y]%輸出
實(shí)驗(yàn)結(jié)果
實(shí)驗(yàn)代碼(僅改動部分)
%--------------------函數(shù)定義-------------------
function z=f1(x,y)
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
z=-50*y+50*x^2+2*x;
end
%--------------------函數(shù)定義-------------------
a=0;b=1;h=0.1; %邊界及步長
x0=a;y0=1/3; %初值
T=[N,x,y] %輸出解
dy=y-(1/3*exp(-50*x)+x.^2) %輸出與精確解的誤差
實(shí)驗(yàn)結(jié)果過多呜达,此處略去
相關(guān)文章:
MATLAB數(shù)值分析之?dāng)?shù)值積分(一)