⑴ 急求MATLAB編程源代碼用四階龍格庫塔法解如下微分方程 y'=y-2x/y(0<x<1),y(0)=1,步長為h=0.2
%以下另存為文件myrk4.m
function[x,y]=myrk4(ufunc,y0,h,a,b)
%參數:函數名稱,初始值向量,步長,時間起點,時間終點
n=floor((b-a)/h);%求步數
x(1)=a;%時間起點
y(:,1)=y0;%賦初值
%按龍格庫塔方法進行求解
forii=1:n
x(ii+1)=x(ii)+h;
k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);
k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);
k4=ufunc(x(ii)+h,y(:,ii)+h*k3);
y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
end
以下是主程序
%y'=y-2x/y(0<x<1),y(0)=1,步長為h=0.2
fun=inline('y-2*x/y');
[t1,f1]=myrk4(fun,1,0.2,0,1);%測試時改變test_fun的函數維數,別忘記改變初始值的維數
subplot(211);plot(t1,f1)%自編函數
title('自編函數求解結果')
%用系統自帶函數ode45進行比較
[t,f]=ode45(fun,[01],1);
subplot(212);plot(t,f);title('ode45求解結果')
⑵ matlab這個代碼怎麼寫
使用matlab計算該遞推數列,需要要符號計算symbolic toolbox以提高計算精度,從而實現數值方法。對於精度,取|pₙ₊₁-pₙ|<ε為迭代終止條件,計算過程為:依次計算aₙ,bₙ,cₙ,sₙ,pₙ,並檢查|pₙ-pₙ₋₁|<ε,一符合即結束循環,代碼如下:
digits(30);
ab0=[sym(1) ,sym(1/2)^0.5 sym(0)];
r=@(x)(x(1:end-1));
n=@(x)(x(end));
abab=@(t)([mean(r(t)) prod(r(t))^0.5 n(t)+1]);
ab1=abab(ab0);
c=@(t)(sum((r(t).^2).*[1 -1]));
s0=sym([1]);
ss=@(t,x)(t-c(x).*2*2^n(x));
s1=ss(s0,ab1);
ee=sym(10e-19);
p=@(ab,s)(sum(r(ab)).^2/4/s);
p0=p(ab0,s0);
p1=p(ab1,s1);
pold=1000;
while abs(pold-p1)>=ee;
fprintf('p(%d)=',n(ab1));
disp(vpa(p1,29));
ab1=abab(ab1);
s1=ss(s1,ab1);
pold=p1;
p1=p(ab1,s1);
end;
結果如上
⑶ 求MATLAB代碼
我只能說,這種限定行數的東西很疼。
再說句不好聽的,多一些空行,代碼一行拆成兩行寫,多一些注釋,再牛一點把一些MATLAB自帶的函數用循環再實現一次,這樣的話,幾百行的代碼肯定可以干到1200行。