超感動，簡單matlab範例，終於做出來了，Oh YA！
經過無數次失敗，在不同的地方犯錯，終於，感動的一刻，不知不覺地來到。
呼，終於可以挑戰paper用的複雜公式了。
Finite difference真的是很可怕，錯一步後，就會給你很好笑的解。
範例題：one-dimensional and time-dependent diffusion equation.
以下是用10.34老師Beers的檔案去改的。解這個問題更簡單的方法，可以直接在matlab下指令edit pdex1，只需要輸入一些參數，不需要自己寫程式，很簡單，但缺乏修改的彈性。
function iflag_main = diffusion_stiffness(N);
clear all; clc; close all;
iflag_main = 0;
% generate 1-D grid
N=40; x = linspace(0,1,N+2); dx = x(2) - x(1);
% set initial condition
for i=1:N
phi0(1,i) = 1*sin(pi*x(i+1));
end
% set vector of times at which to print answer
t_end = 0.5; t_space = 0.1;t_vect = [0:t_space:t_end]
% First simulate using ode15s solver with BDF method
OptionsODE = odeset('BDF','on');
cpu_ode15s_BDF = cputime;
[t,phi] = ode15s(@diff_calc_f,t_vect,phi0,OptionsODE,N);
cpu_ode15s_BDF = cputime - cpu_ode15s_BDF;
make_phi_plot(t,x,phi,'ode15s-BDF',cpu_ode15s_BDF,N, t_space);
% Next, do simulation again with BDF option off
cpu_ode15s_default = cputime;
[t,phi] = ode15s(@diff_calc_f,t_vect,phi0,OptionsODE,N);
cpu_ode15s_default = cputime - cpu_ode15s_default;
make_phi_plot(t,x,phi,'ode15s-default',cpu_ode15s_default,N, t_space );
% compute various sim times compared to ode15s-BDF
disp(' '); disp('Relative simulation times');disp('ode15s-BDF, relative CPU time = 1');
disp(['ode15s-default, relative CPU time = ', ... num2str(cpu_ode15s_default/cpu_ode15s_BDF)]);
u = exp(-t)*sin(pi*x); % analytic solutionplot(x,u,'rx');
iflag_main = 1;
return;
%--------------------------------------------------------
% This MATLAB routine returns the time derivatives of
% each grid point value for the 1-D diffusion problem.
% K.J. Beers. MIT ChE. 9/12/03
function f = diff_calc_f(t,phi,N);
% generate 1-D grid
t; x = linspace(0,1,N+2); dx = x(2) - x(1);
% set A matrix
A = spalloc(N,N,3*N);A(1,1) = -2; A(1,2) = 1;
for k=2:(N-1)
A(k,k-1) = 1; A(k,k) = -2; A(k,k+1) = 1;
end
A(N,N-1) = 2; A(N,N) = -2;
Jac = A./(dx^2)/pi^2; % set Jacobian matrix of system
b=zeros(N,1); b(N,1)=-2*dx*pi*exp(-t)/(dx^2)/pi^2;
f = Jac*phi+b;
return;
%---------------------------------------------------------
% This MATLAB function makes a plot of the phi vs. x
% at various times during the simulation.
% K.J. beers. MTI ChE 9/12/03
function iflag = make_phi_plot(t,x,phi,method_name,cpu_time,N, t_space );
iflag = 0;
x = linspace(0,1,N+2);dx = x(2) - x(1); phi_re=zeros(length(t),N+2);
for k=1:length(t)
for i=1:N
phi_re(k,i+1)=phi(k,i);
end
phi_re(k,N+2)=phi_re(k,N)-2*dx*pi*exp(-t_space*(k-1));
end
figure;for k=1:length(t)
plot(x,phi_re(k,:)); hold on;
end
xlabel('x'); ylabel('\phi(t,x)');
title(['FD solution of time-dependent diffusion problem: ', ... method_name, ', CPU = ', num2str(cpu_time)]);
iflag = 1;
return;
Read more...
Posts
商業週刊的封面：讓我想到I-Fan在某次沙龍裡講藝術竊盜案，提到的吶喊(大家都知道這幅畫吧)，據說有四幅差不多的真畫(是這樣嗎？記性變差了...)，而且重點是，不是畫裡的人在吶喊(雖然我們都這樣覺得)，而是畫裡的人聽到有人在吶喊，覺得很吵，所以摀住自己的耳朵(有沒有人本來就知道的啊？)，總之，很傳神的一幅畫，超好用，像JOBS PANIC(工作恐慌)，到處都在哀嚎，啊~~~
ps. 除了我的advisor的笑容很迷人之外，最近發現ISN的某visiting scientist的笑容也很迷人耶。雖然，本篇過於著重可愛的外國人，但是，我應該是不會跟外國人交往的啦(哈)。
