龍格庫塔解衰減方程的演算法

㈠ matlab 龍格庫塔解的復雜微分方程組，自己編的計算簡單的方程可以，解復雜的方程組可能不對。請高手指點

用matlab的龍格-庫塔法求解復雜微分方程組的關鍵，在於建立自定義微分方程組函數。你的問題主要是沒有利用數組格式來寫，沒有理順其之間的關系。具體應按下列格式書寫

function dy=odefun(t,y)

ft=0.005*t^3-350*t^2+0.002*t;

hy1=0.5*y(1)^3+0.002*y(1)^2+0.005;

dy(1)=500*ft*sin(t)/(hy1*y(1));

dy(2)=y(2)*(y(1)+(700-2*y(2))*tan(0.5*t))/(700-y(1)*tan(0.5*t));

dy=dy(:);

end

調用格式為

y0=[0.001,0];

[t,y]=ode45(@odefun,[0.58:-0.001:0],y0)

運行結果

㈡ matlab龍格庫塔法解微分方程

function [Y] = RK45(t,X,f,h)
K1=f(t,X);
K2=f(t+h/2,X+h/2*K1);
K3=f(t+h/2,X+h/2*K2);
K4=f(t+h,X+h*K3);
Y=X+h/6*(K1+2*K2+2*K3+K4);
end
以上是4階龍格庫塔法的代碼:
自己寫函數，存為f.m
function dxdt = f (t,x)
dxdt(1)=exp(x(1)*sin(t))+x(2);
dxdt(2)=exp(x(2)*cos(t))+x(1); % x(1)是你的f，x(2)是你的g
dxdt=dxdt(:);
end

自己給出t0,x0,h的值（初始時間，初值，步長）
如果求t0到t1的軌跡的話:給個例子如下

t0=0;t1=5;h=0.02;x0=[-1;-1];
T=t0:h:t1;X=zeros(length(x0),length(T));X(:,1)=x0;
for j=1:length(T)-1
X(:,j+1)=RK45(T(j),X(:,j),@(t,x) f(t,x),h);
end
plot(T,X(1,:));
hold on;
plot(T,X(2,:),'r');
具體參數自己設置

㈢ 用 四階龍格庫塔 解三階微分方程 公式和程序

初值給的不夠啊。
不僅要給y,還要給y的一階導和二階導。否則數值解沒法弄。

理論解的話，可能還可以含有一些系數。

==============================================================
在Matlab下輸入：edit，然後將下面兩行百分號之間的內容，復制進去，保存
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DYDt=_amanm(t,Y)
y=Y(1);
dydt=Y(2);
d2ydt2=Y(3);
DYDt=[Y(2);Y(3);(10-10*y-dydt-0.11d2ydt2)/0.001];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

在Matlab命令行下面輸入：
t_start=0;
t_end=1;
yy0=[0.1;0;0]; %初值
[t,y]=ode45('_amanm',[0,t_end],yy0);
plot(t,y(:,1));
xlabel('t');
ylabel('y')
y(end,1)
得到的結果：
ans =

1.0056

用昨天給的C代碼計算出來的結果是：
1.005630

㈣ mathematic用龍格庫塔法解二階方程

……所謂龍格庫塔法，通俗地說，就是把一個n階的常微分方程，整理成n個形如 f'(t)=g(t,f(t)) （注意此時右側不含 f(t) 的導數）的一階常微分方程組再加以求解的方法。你的方程整理成龍格庫塔所需要的形式就是：

x'[t] = y[t]

y'[t] = (-c y[t] - k x[t] - F[t])/m

再考慮到你找到的這段代碼本身對格式的要求，就有：

m=1;c=1;k=1;f=1;
stepsize=0.01;t=1;

sol=NestList[
RungeKutta[{y,(-cy-kx-f)/m},#,{x,y},stepsize]&,{0.,0.},Round[t/stepsize]];

剩下的你自己弄吧……算了，考慮到你的方程里的F[t]涉及了t的具體值，又多了一點點坎（在要使用你找到的這個函數的前提下），再來一個F[t_]=Sin[t]的代碼示例吧：

Clear[f,t]
m=1;c=1;k=1;
f[t_]=Sin[t];
stepsize=0.01;tmin=0;tmax=1;

sol=FoldList[
RungeKutta[{y,(-cy-kx-f[#2])/m},#,{x,y},
stepsize]&,{0.,0.},Range[tmin,tmax,stepsize]];

㈤ 龍格庫塔的文字介紹

龍格-庫塔（Runge-Kutta）法
到目前為止，我們已經學習了多步法，例如：亞當斯-巴什福思（Adams
-Bashorth）法，亞當斯-莫爾頓（Adams-Monlton）法，都是常微分方程的積分
方法。它們需要在每一次迭代時重新計算一遍等式右邊的結果（非線性隱含問
題忽略計算多個 f (ω)值的可能性）
龍格-庫塔（Runge-Kutta）法是一種不同的處理，作為多級方法為人們所知。
它要求對於一個簡單的校正計算多個 f 的值。
下面，我們列出了 3 種最流行的龍格-庫塔（Runge-Kutta）法：
改進的歐拉方法（精度：p=2）：
V a = V n + Δtf (V n,tn)
2
Δt)二階格式
V n+1 = V n +Δtf (V a,tn +
2
Hevn』s 方法（p=2）：
這是另一種二階格式：
V a = V n +Δtf (V n,tn)
V n = V n +
+1 Δt[ f (V n,tn) + f (V a,tn +Δt)]
2
注意： f (Vn,tn)在運算中應該只被計算一次。
四次龍格-庫塔（Runge-Kutta）法（p=4）：
這是一個 4 階格式。這次我們寫的形式有點不同：
a = Δtf (V n,tn)
b = Δtf (V n + 1 a,tn + 1
2 2 Δt)
c = Δtf (V n + 1 b,tn + Δt)
1
2 2
d = Δtf (V n + c,tn +Δt)
V n =V n +
+1 1 (a +2b +2c +d)。
6

㈥ matlab龍格庫塔法求解微分方程，怎麼編程

用matlab龍格庫塔法（ode45）求解微分方程組的方法。

第一步：建立微分方程阻的自定義函數myodefun（t，z），其主要內容

dz(1)=1.2*z(1)-0.3*z(1)*z(2);

dz(2)=0.05*z(1)*z(2)-0.1*z(2);

第二步：用ode45函數求解。

[t,z]=ode45(@myodefun,[0,30],[10,3])

第三步：用plot函數繪制t-z1(t)和t-z2(t)的圖形。

計算結果

㈦ 龍格庫塔求解二階微分方程組的MATLAB編程

MATLAB求解x''+0.7x'+0.8x'|x'|+25.6x-25.6x³=0二階微分方程組的方法，可以按下列步驟進行：
1、建立自定義函數func（）
function
f
=
func(t,x)
%x''+0.7x'+0.8x'|x'|+25.6x-25.6x³=0
f(1)=x(2);
f(2)=25.6*x(1)^3-25.6*x(1)-0.8*x(2)*abs(x(2))-0.7*x(2);
f=f(:);
2、建立龍格庫塔演算法函數runge_kutta（）
調用格式：[t,x]
=
runge_kutta(@(t,x)func(t,x),x0,h,a,b);
3、然後根據x和x'數據，繪制出x(t)、x′(t)的圖形。
plot（x（：，1）,x（：，2））

㈧ MATLAB中已知系統微分方程及初始值用歐拉法和龍格庫塔法解一階微分方程

function Euler
%歐拉法和龍格庫塔演算法解一階常微分方程源代碼
%例子dy/dx=-y+x+1
f=inline('-y+x+1','x','y'); %微分方程的右邊項
dx=0.5; %x方向步長
xleft=0; %區域的左邊界
xright=10; %區域的右邊界
xx=xleft:dx:xright; %一系列離散的點
n=length(xx); %點的個數
y0=1;
%%(1)歐拉法
Euler=y0;
for i=2:n
Euler(i)=Euler(i-1)+dx*f(xx(i-1),Euler(i-1));
end
%%(2)龍格庫塔法
RK=y0;
for i=2:n
k1=f(xx(i-1),RK(i-1));
k2=f(xx(i-1)+dx/2,RK(i-1)+k1*dx/2);
k3=f(xx(i-1)+dx/2,RK(i-1)+k2*dx/2);
k4=f(xx(i-1)+dx,RK(i-1)+k3*dx);
RK(i)=RK(i-1)+dx*(k1+2*k2+2*k3+k4)/6;
end
%%Euler和Rk法結果比較
plot(xx,Euler,xx,RK)
hold on
%精確解用作圖
syms x
rightsolve=dsolve('Dy=-y+x+1','y(0)=1','x');%求出解析解
rightdata=subs(rightsolve,xx);%將xx代入解析解,得到解析解對應的數值
plot(xx,rightdata,'r*')
legend('Euler','Runge-Kutta','analytic')

㈨ 龍格庫塔求解微分方程組

只要理解了龍格庫塔，這就很容易了。
定義函數，
f=func(x,y)
if y(1)>126,
f=[f1(y(1),y(2)), f2(y(1),y(2))];
else
f=[g1(y(1),y(2)), g2(y(1),y(2))];
end

然後就簡單了

㈩ 龍格庫塔法的基本原理

該演算法是構建在數學支持的基礎之上的。對於一階精度的拉格朗日中值定理有：
對於微分方程：y'=f(x,y)
y(i+1)=y(i)+h*K1
K1=f(xi,yi)
當用點xi處的斜率近似值K1與右端點xi+1處的斜率K2的算術平均值作為平均斜率K*的近似值，那麼就會得到二階精度的改進拉格朗日中值定理：
y(i+1)=y(i)+[h*( K1+ K2)/2]
K1=f(xi,yi)
K2=f(x(i)+h,y(i)+h*K1)
依次類推，如果在區間[xi,xi+1]內多預估幾個點上的斜率值K1、K2、……Km，並用他們的加權平均數作為平均斜率K*的近似值，顯然能構造出具有很高精度的高階計算公式。經數學推導、求解，可以得出四階龍格－庫塔公式，也就是在工程中應用廣泛的經典龍格－庫塔演算法：
y(i+1)=y(i)+h*( K1+ 2*K2 +2*K3+ K4)/6
K1=f(x(i),y(i))
K2=f(x(i)+h/2,y(i)+h*K1/2)
K3=f(x(i)+h/2,y(i)+h*K2/2)
K4=f(x(i)+h,y(i)+h*K3)
通常所說的龍格-庫塔法是指四階而言的，我們可以仿二階、三階的情形推導出常用的標准四階龍格-庫塔法公式