微分方程組的四階龍格庫(kù)塔公式求解matlab版
微分方程組的四階龍格庫(kù)塔公式，求解matlab版的代碼。下面是我的代碼，運(yùn)行不出來(lái)。
function varargout=rungekutta(varargin)
clc,clear
x0=0;xn=1.0;y0=1;h=0.05;%h為步長(zhǎng)
[y,x]=rungekutta4(x0,xn,y0,h);
Function z=f(x,y);
z=y-2*x/y;
function [y,x]=rungekutta4(x0,xn,y0,h)
x=x0:h:xn;
n=(xn-x0)/h;
y1=x;
y1(1)=y0;
for i=1:n %龍格庫(kù)塔法的算法
K1=f(x(i),y1(i));
K2=f(x(i)+h/2,y1(i)+h/2*K1);
K3= f(x(i)+h/2,y1(i)+h/2*K2);
K4= f(x(i)+h,y1(i)+h*K3);
y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);
yy1(i)=(1+2*x(i+1)) ^0.5;
error(i)=y1(i+1)-yy1(i);
yy(i+1)=yy1(i);
e(i+1)=error(i);
end
y=y1;
n=(xn-x0)/h;
fprintf(‘i x(i) yi數(shù)值解 y(i)真值 誤差\n’);
fprintf(‘-------------------------------------------------------------\n’);
for i=1:n
fprintf(‘%2d %12.4f %14.8f %14.8f %12.8f\n’,i, x(i+1), y(i+1), yy(i+1), error(i));
end
plot(x,e,‘r.’);%畫(huà)出誤差曲線(xiàn)
xlable(‘x軸’),ylable(‘誤差’);
title(‘步長(zhǎng)為0.05時(shí)的誤差曲線(xiàn)’);
返回小木蟲(chóng)查看更多
京公網(wǎng)安備 11010802022153號(hào)
一、幾個(gè)自定義函數(shù)位置錯(cuò)誤。應(yīng)為這樣放置
1、主程序 ，
function rungekutta()
clc,clear
x0=0;xn=1.0;y0=1;h=0.05;%h為步長(zhǎng)
[y,x]=rungekutta4(x0,xn,y0,h);
end
2、 四階龍格庫(kù)塔函數(shù)程序
function [y,x]=rungekutta4(x0,xn,y0,h)
。。。
end
3、z函數(shù)程序
function z=f(x,y);
z=y-2*x/y;
end
二、xlable(‘x軸’),ylable(‘誤差’);這句代碼中函數(shù)書(shū)寫(xiě)錯(cuò)誤。應(yīng)為
xlabel('x軸'),ylabel('誤差');
三、修改后運(yùn)行得到如下結(jié)果
什么結(jié)果？大神。
運(yùn)行結(jié)果：
i x(i) yi數(shù)值解 y(i)真值 誤差
-------------------------------------------------------------
1 0.0500 1.04880886 1.04880885 0.00000001
2 0.1000 1.09544514 1.09544512 0.00000002
3 0.1500 1.14017546 1.14017543 0.00000004
4 0.2000 1.18321600 1.18321596 0.00000005
5 0.2500 1.22474493 1.22474487 0.00000006
6 0.3000 1.26491113 1.26491106 0.00000007
7 0.3500 1.30384056 1.30384048 0.00000008
8 0.4000 1.34164088 1.34164079 0.00000009
9 0.4500 1.37840498 1.37840488 0.00000011
10 0.5000 1.41421368 1.41421356 0.00000012
11 0.5500 1.44913781 1.44913767 0.00000014
12 0.6000 1.48323985 1.48323970 0.00000015
13 0.6500 1.51657526 1.51657509 0.00000017
14 0.7000 1.54919353 1.54919334 0.00000019
15 0.7500 1.58113904 1.58113883 0.00000021
16 0.8000 1.61245178 1.61245155 0.00000023
17 0.8500 1.64316793 1.64316767 0.00000026
18 0.9000 1.67332034 1.67332005 0.00000028
19 0.9500 1.70293895 1.70293864 0.00000031
20 1.0000 1.73205115 1.73205081 0.00000034
圖像上傳不上來(lái)
，
function rungekutta()
clc,clear
x0=0;xn=1.0;y0=1;h=0.05;%h為步長(zhǎng)
[y,x]=rungekutta4(x0,xn,y0,h);
End
function [y,x]=rungekutta4(x0,xn,y0,h)
x=x0:h:xn;
n=(xn-x0)/h;
y1=x;
y1(1)=y0;
for i=1:n %龍格庫(kù)塔法的算法
K1=f(x(i),y1(i));
K2=f(x(i)+h/2,y1(i)+h/2*K1);
K3= f(x(i)+h/2,y1(i)+h/2*K2);
K4= f(x(i)+h,y1(i)+h*K3);
y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);
yy1(i)=(1+2*x(i+1))^0.5;
error(i)=y1(i+1)-yy1(i);
yy(i+1)=yy1(i);
e(i+1)=error(i);
end
function z=f(x,y);
z=y-2*x/y;
end
y=y1;
n=(xn-x0)/h;
fprintf('i x(i) yi數(shù)值解 y(i)真值 誤差\n');
fprintf('-------------------------------------------------------------\n');
for i=1:n
fprintf(‘%2d %12.4f %14.8f %14.8f %12.8f\n’,i, x(i+1), y(i+1), yy(i+1), error(i));
end
plot(x,e, 'r. ');%畫(huà)出誤差曲線(xiàn)
xlabel('x軸'),ylabel('誤差');
title('步長(zhǎng)為0.05時(shí)的誤差曲線(xiàn)');
大神，我修改完，還是運(yùn)行不出來(lái)啊。你能不能直接把你的完整代碼發(fā)給我，我自己找找錯(cuò)誤。
怎么發(fā)給你？