Matlab 非線性熱傳導(拋物方程)問題

2022-08-03 20:36:20 字數 3230 閱讀 5678

函式檔案1:real_fun.m

1 function f=real_fun(x0,t0)

2 %精確解

3 f=4*x0*(1-x0)*sin(t0);

函式檔案2:f.m

1 function f=f(n,u,u,t,h1,h2)

2 %非線性方程組

3 %h1是x的步長,h2是t的步長

4 %u表示迭代節點,上一時刻的數值解

5 %h表示時間節點上的步長

6 %n表示空間節點的步數

7 a0=0.5*t^4*h2*n^2;

8 f(1,1)=a0*(u(2)^2-2*u(1)^2)+h2*fi(h1,t)+u(1)-u(1);

9 f(n-1,1)=a0*(-2*u(n-1)^2+u(n-2)^2)+h2*fi((n-1)*h1,t)+u(n-1)-u(n-1);

10for p=2:n-2

11 f(p,1)=a0*(u(p+1)^2-2*u(p)^2+u(p-1)^2)+h2*fi(p*h1,t)+u(p)-u(p);

12 end

函式檔案3:fi.m

1 function f=fi(x0,t0)

2 %等式右邊的f函式

3 f=4*x0*(1-x0)*cos(t0)-16*t0^4*(6*x0^2-6*x0+1)*(sin(t0))^2;

函式檔案4:jacobian.m

1 function g=jacobian(n,u,t,h1,h2)

2 %計算每個時間節點的牛頓迭代過程中的雅可比矩陣

3 %u表示迭代初值,上一時刻的數值解作為迭代初值

4 a=0.5*t^4*h2*n^2;

5 g=zeros(n-1);

6 g(1,2)=2*a*u(2);

7 g(1,1)=-4*a*u(1);

8 g(n-1,n-1)=-4*a*u(n-1);

9 g(n-1,n-2)=2*a*u(n-2);

10for p=2:n-2

11 g(p,p+1)=2*a*u(p+1);

12 g(p,p)=-4*a*u(p);

13 g(p,p-1)=2*a*u(p-1);

14end

15 g=g-eye(n-1);

函式檔案5:newtond.m

1 function x=newtond(n,u,t,h1,h2)

2 %使用修改後的牛頓迭代,可以不求雅可比de逆

3 %u中間代初值

4 %u起始迭代初值

5 u=u;

6 tol=0.5e-5;

7 % jacobi=jacobian(n,u,t,h1,h2);%每隔k步求一次雅可比

8 x1=u-jacobian(n,u,t,h1,h2)\f(n,u,u,t,h1,h2);

9while (norm(x1-u,1)>=tol)

10 %數值解的1範數是否在誤差範圍內

11 u=x1;

12 x1=u-jacobian(n,u,t,h1,h2)\f(n,u,u,t,h1,h2);

13end

14 x=x1;%不動點

1

tic;

2clc

3clear

4 n=100;

5 m=1000;

6 t_h=1/m;%t的步長

7 x_h=1/n;%x的步長

8 x=0:x_h:1;%x的節點

9 ti=0:t_h:0.5;%t的節點

10 %********************真解**************************

11for i=1:length(x)

12for j=1:length(ti)

13 real_z(i,j)=real_fun(x(i),ti(j));

14end

15end

16 %********************真解**************************

17 %********************數值解**************************

18 ui=zeros(length(x)-2,1);%牛頓迭代初值

19 z=zeros(length(x),length(ti));

20for i=1:length(ti)-1

21 z(2:length(x)-1,i+1)=newtond(length(x)-1,ui,ti(i+1),x_h,t_h);%t(i+1)時間的牛頓數值解

22 ui=z(2:length(x)-1,i+1);%牛頓迭代初值,上一時刻的數值解作為迭代初值

23end

2425 %********************數值解**************************

26 [x,y]=meshgrid(x,ti);

27 subplot(2,2,1),

28 mesh(x,y,real_z');

29 xlabel('x');ylabel('t');zlabel('u');title('analytical solution');

30 subplot(2,2,2),

31 mesh(x,y,z');

32 xlabel('x');ylabel('t');zlabel('u');title('numerical solution');

33 subplot(2,2,3),

34 mesh(x,y,real_z'-z');

35 xlabel('x');ylabel('t');zlabel('u');title('error solution');

36 title('牛頓迭代法');

37grid on;

38 toc;

效果圖:

求解熱傳導方程matlab

這是非穩態一維熱傳導的方法,也叫古典顯格式。如果是做數學建模,就別用了,這種方法計算量比較大,算的很慢,而且收斂不好。但是如果實在沒辦法也能湊合用。該改的地方我都用?代替了。給個詳細解釋 1 function rechuandao 23 llist 4 n 1000 空間點數 5 m 10000 6...

matlab 非線性優化

求解非線性問題 min z f x s.t.c x 0,ceqx 0,ax b,aeqx beq,lb x ub.x,fval,exitflag,output,lambda,grad,hessian fmincon fun,x0,a,b,aeq,beq,lb,ub,nonlcon,options,p...

Matlab非線性規劃

在matlab非線性規劃數學模型可以寫成一下形式 minf x s.t.begin ax le b aeq x beq c x le 0 ceq x 0 end f x 為目標函式,a,b,aeq,beq為線性約束對應的矩陣和向量,c x ceq x 為非線性約束。matlab求解命令為 x fmi...