簡潔版
code
fprintf('請輸入區間下界:\n')
a=input('');
fprintf('請輸入區間上界:\n')
b= input('');
fprintf('請輸入初值alpha:\n')
alpha=input('');
fprintf('請輸入最大迭代次數n:\n')
n = input('');
x0=a;
y0=alpha;
h=(b-a)/n;
k=;x1=;
y1=;
for n=1:n
k1=h*fun_kutta(x0,y0);
k2=h*fun_kutta(x0+h/2,y0+k1/2);
k3=h*fun_kutta(x0+h/2,y0+k2/2);
k4=h*fun_kutta(x0+h,y0+k3/2);
x1(n)=x0+h;
y1(n)=y0+(k1+2*k2+2*k3+k4)/6;
x0=x1(n);
y0=y1(n);
end
function z=fun_kutta(x,y)
% z = x + y;
% z = - y^2;
% z=2*y/x+x^2*exp(x);
% z=(y^2+y)/x;
% z=-20*(y-x^2)+2*x;
% z=-20*y+20*sin(x)+cos(x);
z=-20*(y-exp(x)*sin(x))+exp(x)*(sin(x)+cos(x));
end
fprintf('請輸入區間下界:\n')
a=input('');
fprintf('請輸入區間上界:\n')
b= input('');
fprintf('請輸入初值alpha:\n')
alpha=input('');
% fprintf('請輸入最大迭代次數n:\n')
% n = input('');
k=;x1=;
y1=;
fp=fopen('3.3.1.txt','a');%'a.txt'為檔名；'a'為開啟方式：在開啟的檔案末端新增資料，若檔案不存在則建立。
for n=5:5:20
x0=a;
y0=alpha;
h=(b-a)/n;
for n=1:n
k1=h*fun_kutta(x0,y0);
k2=h*fun_kutta(x0+h/2,y0+k1/2);
k3=h*fun_kutta(x0+h/2,y0+k2/2);
k4=h*fun_kutta(x0+h,y0+k3/2);
x1(n)=x0+h;
y1(n)=y0+(k1+2*k2+2*k3+k4)/6;
x0=x1(n);
y0=y1(n);
endfprintf(fp,'\n');
fprintf(fp,'%f\t',x1);
fprintf(fp,'\n');
fprintf(fp,'%f\t',y1);%fp為檔案控制代碼，指定要寫入資料的檔案。注意：%d後有空格。
fprintf(fp,'\n');
end fclose(fp);%關閉檔案。
四階龍格庫塔法
這裡主要講一下如何用c語言程式設計運用四階龍格庫塔法求解微分方程組。對於所舉例子，只是為了說明龍格庫塔法不僅可以解一階線性微分方程，高階非線性也可通過降階後按照經典四階龍格庫塔法公式逐步求解。只要選取合適的步長h，就能夠平衡速度和精度，達到求解要求。至於例子中的一級倒立擺的物理含義沒有提及到，各種方...
四階龍格庫塔法的基本思想 四階龍格庫塔實驗報告
1 三 四階runge kutta法求解常微分方程 一 龍格庫塔法的思想根據第九章的知識可知道，euler方法的區域性截斷誤差是，而當用euler方法估計出再用梯形公式進行校正，即採用改進euler方法得出數值解的截斷誤差為。由lagrange微分中值定理記，得到這樣只要給出一種計算的演算法，就能得...
python實現四階龍格庫塔法
coding utf 8 created on sun dec 24 15 29 08 2017 author www 本程式是用四階龍格庫塔法求解課本 數值計算方法 馬東昇 p242頁的例7 3 fun為指定的導數的函式 rf4為四階龍格庫塔法 def fun x,y f y 2 x y retu...