Here is the problem
Ive got a FEM code
function main1clc n = 100; x0 = ones(3*n,1); sol = fsolve(@(x)fun(x,n),x0); norm(fun(sol,n)) x = ((1:n)-1)/(n-1); plot(x,sol(1:n))endfunction res = fun(z,n) eta=1.0; beta = 1.0; x = ((1:n)-1)/(n-1); dx = 1/(n-1); y1 = z(1:n); y2 = z(n+1:2*n); y3 = z(2*n+1:3*n); for i=1:length(y1) Y1=y1(i); end for i=1:length(y2) Y2=y2(i); end res_y1 = zeros(n,1); res_y2 = zeros(n,1); res_y3 = zeros(n,1); res_y1(1) = y1(1)-10.0; for i=2:n-1 res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/Y2); end res_y1(n) = y1(n); res_y2(1) = y2(1); for i = 2:n-1 res_y2(i) = (y2(i+1)-2*y2(i)+y2(i-1))/dx^2 - y1(i).^2; end res_y2(n) = y2(n)-eta; res_y3(1) = y3(1); for i=2:n-1 res_y3(i) = (y3(i+1)-2*y3(i)+y3(i-1))/dx^2 - Y2; end res_y3(n) = y3(n)-1.0; res = [res_y1;res_y2;res_y3];end
It seems that using
res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/Y2);
and
res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/y2(i));
must be eqiuvalent, but NO. The results differ
Best Answer